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THE INTERACTION BETWEEN A SOLID BODY AND VISCOUS 
FLUID BY .MARKER-AND-CELL METHOD 

By 

Rob er t Y.K. Cheng 1 
SUMMARY 


A computational method for solving nonlinear problems relating 
to impact and penetration of a rigid body into a fluid-type medium 
is presented. The numerical technique, based on the Marker-and-Cell 
method, gives the pressure and velocity of the flow field. An impor- 
tant feature in this method is that the force and displacement of 
the rigid body interacting with the fluid during the impact and 
sinking phases are evaluated from the boundary stresses imposed 
by the fluid on the rigid body. 

A sample problem of low velocity penetration of a rigid block 
into still water is solved by this method. The computed time his- 
tories of the acceleration, pressure, and displacement of the block 
show good agreement with experimental measurements. A sample problem 
of high velocity impact of a rigid block into soft clay is also 
presented. The constitutive relationships of the clay is repre- 
sented as a very viscous non-Newtonian fluid. 


1 Professor of Civil Engineering, School of Engineering, Old 
Dominion University, Norfolk, Virginia 23508. 



INTRODUCTION 


There is a need to develop an analytical method of studying 
the problem of aircraft landing impact on soil runways. Since 
deep ruts are created by the tires on unpaved runways, the drag 
forces on the landing ( gear can be large enough to endanger the 
safe operation of the aircraft. 

Many studies have been performed {refs. 1, 2, 3, and 4) to 
study the soil-tire interactions at ground speeds associated with 
aircraft operations. The methods used in those studies are based 
on either empirical modeling techniques or on soil strengths 
obtained by static tests. Although the rut depth and drag forces 
will depend on many factors other than the soil strength, a rational 
method for predicting the landing- impact problem must, take into 
account the dynamic soil properties, the interface- stress distri- 
bution between the tire and soil, and .the large movement of soil 
masses. 

As a first step in a program to provide an analytic method 
as an engineering tool for studying the landing- impact problem, 
a computational method is presented for solving the problem of 
large displacements of fluid-acting media when interacting with 
a solid body. An important feature in this method is that the 
force and displacement of the rigid body during the impact and 
sinking phases are derived from the fluid field. 

LIST OF SYMBOLS 

C wave speed 

D divergence of fluid medium 

E energy 

e strain-rate 

e strain-rate tensor 

f force component contributed by a single cell • 
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force component acting on rigid body 
force vector 

gravitational acceleration vector 
gravitational acceleration 

depth of fluid field measured from the base of rigid body 

unity dyadic 

length of fluid field 

mass of rigid body 

direction cosine of unit tangential vector 
direction cosine of unit normal vector 
hydrostatic pressure 
velocity vector 

source term of Poisson's equation for calculating <f> 
time 

velocity component in x-direction 

velocity component in y-direction 

cartesian coordinates (also used as subscripts) 
time increment 

width and height, respectively, of all cells 
strain 

strain tensor 
normal stress 
stress tensor 
shear stress 

deviatoric stress tensor 

del operator 

dilation 
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y modulus of rigidity 

tc bulk modulus 

r) non-Newtonian viscosity 

apparent viscosity 
v kinematic viscosity 

<p ratio of pressure to mass density of fluid 

p mass density of fluid 

X convergence criterion 

X f force tolerance for iterative scheme 


Subscripts 


11 , 22/22 

e 


i/j 

o 

w 

a 


diagonal elements of stress tensor 
identifies trial value for rigid body 
denotes x and y direction, respectively 
initial impact velocity 
identifies rigid body 
externally applied pressure 


Superscripts 

" spherical stress tensor 

5 deviatoric stress tensor 

n time cycle 

m order of recycle 

s order of iteration 



SOIL PROPERTIES AND CONSTITUTIVE RELATIONSHIPS 


A valid solution in continuum mechanics requires that the 
solution satisfies the equations of conservation of mass and of 
energy, equations of motion, and the equation of state of the 
material. When the tire of an aircraft landing at relatively 

l 

high speed strikes the free surface of the soil, a strain- rate 
which depends upon the landing velocity' is imposed on the wheel 
and the soil. Assuming the deformation will be large, the momentum 
(or energy) of the aircraft will be dissipated by the soil in the 
form of elastic deformations and plastic flow. The elastic defor- 
mations can be expressed as a function of elastic stresses whereas 
the plastic flow can be expressed as functions of shear (or viscous) 
stresses with magnitudes depending on the rate of strain. The 
yield stress will be defined as the transition point from the 
elastic to the plastic state. Plastic flow will appear only 
when the stresses exceed a certain limit indicated by the yield 
stress. For -materials exhibiting distinct yield points, there 
is no ambiguity in establishing this limit. The stress-strain 
behavior of most soils seldom indicate any distinct point that 
may be identified as a yield point. Therefore, in this paper 
the yield stress will be defined as the transition point from the 
linear stress-strain behavior to the hon-linear stress-strain 
behavior. Below the yield point, the stress-strain relationship 
is not rate dependent. 

The various phenomena which are generated by the landing 
impact of an aircraft may be treated as follows: 

a. In the extremely "close-in region" in the neighborhood 

of the impact in which the stresses are very large, with resulting 
large displacements and velocities of the material, the medium 
will be in a state of plastic flow. The impact momentum is pri- 
marily dissipated by' the shear resistance of the material. The 
phenomenon is of a deviatoric nature. 

b. As the magnitude of stresses decreases from the zone of 
loading, there is a transition region in which the medium undergoes 
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a transition from the plastic flow state to the elastic state. 

The momentum dissipation in this region is relatively small in 
comparison with the flow state. The stresses and displacements 
can be evaluated by the elastic constants of the "solid" material. 

c. As the moving load moves away from the "close-in region”, 
the applied load is released, and the flow state reverts to an 
elastic state. The elastic rebound is a rather common phenomenon 
observed upon unloading in many simple compression (or tension) 
tests. It must be noted that the material, after being subjected 
to plastic flow, reverts to a solid but with elastic and flow 
properties which may be drastically different from its virgin 
state. 

Some of the physical variables of sandy and clayey soils 
influencing mass behavior are: air and water content, grain size 


distribution, and density. The environmental variables affecting 
the mass behavior are confining pressures, stress history, rate of 
loading, and duration of loading. The elastic and flow properties 
of soils are all influenced by these variables. The state in which 
the soil remains elastic is highly dependent upon the confining 
pressure. 

The state of 
tensor o , as: 

stress 

at a 

point will be described by a stress 
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The strain-rate tensor e is related to strain tensor as: 



It is to be noted that the above relationship is valid only where 
the strains are small, so that higher order products of strains 
and rates of strains may be neglected. For a small incremental 
stress-strain relationship, as will be the case considered in 
this formulation, the relationship will be valid. 

The velocity of a material particle at some point and at- time 
t is given by the vector q which is completely described by 
time rate of change of the displacement functions. The strain 
rate tensor defined in terms of displacement functions is? 

e = j (.Vq + qV) 

The total deformation in a soil mass may be decomposed into (a) 
volumetric deformation caused by hydrostatic (or spherical) stress 
components, and (b) deviatoric deformation caused by deviatoric 
stress components. The stress, stain, and strain-rate tensors 
may be decomposed into hydrostatic and deviatoric parts given as: 

a" + ff’ 

e" + i 5 (1) 

e" + e’ 

The deviatoric tensor has the property that the first invariant 
is zero. Generally the stress-strain behavior of a medium for 
the hydrostatic stress field is linear. 

From small incremental stress-strain theory, the volumetric 
strain is the dilation A given by: 

A = V ° q 
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The hydrostatic stress tensor is the mean of the normal stress 
components given by; 

-p 0 0 

0 -p 0 
0 0 -p 

where 

X 

-P = 3 (0 n + CT 2 2 + ff 33) 

If the material is isotropic with a bulk modulus k (not necessarily 
a constant) , a constant rigidity modulus y in the elastic state, 
and a flow parameter n in the plastic flow state, the constitutive 
equations are: 

p = kA 

o' = 2 ye* (elastic region only) 

= 2rie 11 (plastic region only) 

In comparison with plastic flow, the shear deformation is rela- 
tively small in the elastic region. The material behavior in the 
elastic region may be associated with the conservative (recoverable) 
part whereas in the plastic region it is associated with the dissi- 
pative (non-recoverable) part, in which the flow takes place at 
constant volume. In large displacement problems, the deformation 
due to the conservative part will be small in comparison with the 
dissipative part. Henceforth, conservative-part deformation will 
be tentatively neglected. For a cohesive soil such as clay, the 
flow parameter n is non-Newtonian and it is a function of the 
confining pressure and strain-rate. A convenient method to express 
the constitutive relationship for a non-Newtonian material is the 
use of an apparent viscosity, ti , (ref. 5) as shown in figure 
1. The constitutive relationships for the Mississippi Buckshot 
clay of various water contents is given in reference 6. 
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GOVERNING EQUATIONS 


The governing equations describing the motion of a continuous 
medium are: 


Conservation of mass 
Conservation of motion 

(3) 

p— = a:i-V «btc, Conservation of energy 
at h 

E = E(p,p) Equation of state 

where b and c, are the heat flux vector and radiative heat term 

h 

respectively. 

For compressible flow, the equation of state is necessary to 
complete the system of equations for solving the unknowns q, p, 
p, and E. For incompressible flow, the governing equations are 
reduced ±o: 


+ pv • g = o 
Pgf = pg + V • o 


V • q = 0 

P {jt + ^ * V )^ = P? + V • CT 


(4) 


Imposing the stress strain relationship expressed in terms of 
displacement functions, the set of governing equations may be 
solved and the solution is given by these displacement functions 
in time derivative f o' rut. 

For reasons given previously, large deformation is attributed 
to plastic flow so that volumetric deformation due to hydrostatic 
stresses will be neglected. The soil system is considered to flow 

i 

under constant volume under the conditions of soil-tire inter- 
action due to a moving wheel of an aircraft. 
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In .the plastic flow state, the constitutive equation of the 
soil system is: 

a = a” + x* = -pi +' 2rie’ 

The deviatoric strain-rate tensor expressed in terms of displacement 
functions is : 

e* = e - e" = ~ (Vq + qV) - y (V » q)I 


The constitutive equation becomes: 


a =-pI + n (Vq + qV) - y T)(V • q) I 


(5) 

: s 

t j 

Substituting the constitutive equation into Eq„ (4), the equation* 
of motion becomes s 


p (jF + 9 • 79 “ P9 - Vp + ri 5 v (V • q) + V 2 q 




'] 


( 6 ) 


+ (Vn) • (Vq + qV) - y (V • q) (Vti) 


Imposing the condition of incompressibility, the governing 
equations describing the motion of a medium are: 


V • q = 0 

p(^r + q 0 v^g - pg - Vp + nv 2 g + (Vn) 0 (Vq + qV) 


(7) 


where n = n(p/e)® For a special case in which the flow parameter 
D is treated as a constant, the equation of motion in Eq. (7) is 

the well known Navier-Stokes equation. 

< 

In view of the non-linear form of Eq. (7) , the difficulties 
involved in obtaining an analytic form of solution satisfying Eq. 
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(7) and the imposed boundary conditions is formidable. By treating 
the soil as an incompressible and viscous material, the governing 
equations are similar to those in fluid dynamic problems. 

The availability of large, high-speed computers and advanced 
numerical techniques (ref. 7) , provides a powerful tool for solving 
complex non-linear boundary value problems in fluid dynamics. A 
computational method based on the Marker and Cell (MAC) technique 

f 

(ref. 8) is used since the primitive variables of velocity com- 
ponents and pressure are solved directly, and the primitive 
variables are required to relate the interaction between the 
solid body and the fluid medium. 

The Marker and Cell (MAC) method is a computational method 
with visual display for solving problems of time-dependent motions 
of a viscous, incompressible fluid with a, free surface.; Some recent 
workers using 'the MAC method are Donovan (ref. 9) and Viecelli 
(ref. 10). The MAC method requires that the external wall shapes 
be confined to the fixed rectangular cells of the Euler ian mesh. 

In this paper, the MAC method has been adapted and modified to 
handle the fluid-solid interaction problem * involving the moving 
wall of the solid body. The restriction of the stationary wall 
boundary has been overcome by an iterative scheme involving the 
impulse-momentum principle and the translation of coordinates. 

t * 

For a complete detailed description and discussion of the 
MAC method, the reader is urged to consult reference 8. Basically 
the method is a numerical technique for solving problems in incom- 
pressible hydrodynamics containing free surfaces. Using two 
spatial dimensions, the primary dependent variables .are the 
pressure and the two velocity components. The variables are 
computed in each time step and the solution is advanced in small 
time increments. Using the Lagrangian approach to represent the 
fluid motion with time, massless "marker 1 ’ particles are moved in 
each time step to new positions according to the velocity com- 
ponents in their vicinities. The marker particles serve to define 
the free surfaces and they do not enter into the solution of the 
Navier-Stokes equations. 
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As with most other numerical methods which work in small 
incremental form, the MAC method works with a small time cycle. 

The results of the field variables in each cycle act as initial 
conditions of the next time cycle, and the calculations are con- 
tinued for as many cycles as the investigator wishes. 

During the ^mall, but finite, time cycle At, the flow 
parameter n may be treated as a constant within each cell, 
such that Vri = 0. Introducing 

4 = 2. and v - — (8) 

v P P 

Eq. (7) describing the motions of an incompressible fluid in two- 
dimensional Cartesian coordinates are: 


D 


3u , 9v 
3x + ‘ 3y 


0 


(9) 


3u 

3t 3x 


u 



3v , 3 

at + ax 


uv + 


3y 


- i£ 

> ax 


+ V 



+ 


H + 
ay 



a^u \ 

a 

8y 2 / 


( 10 ) 


(ii) 


The kinematic viscosity is considered to be a constant, but the 
numerical method can easily account for variable viscosity which 
may be treated as a function of pressure and velocity components. 

a 3 

Operating on Eq. (10) with and Eq. (11) with und 

allowing the interchange of space and time differentials, one 
obtains the Poisson equation: 


V 2 <f> = = -r 

3x 2 Sy 2 


(12) 
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where 


3 2 U 2 . „ 3 UV , 3 2 v 2 , 3D 3 2 D , 3 2 D 

R - + 2 - A- ■ a ' " ' + — + -5+r - v + 

3 x? 9x9 y 3y 2 9t 3x 2 3y 2 

Equations (9) through (12) constitute the basic equations from 
which a finite difference scheme is developed. 

The boundary conditions required fpr the problem are indicated 
in sketch 1. The stationary boundaries are represented by no-slip 



Moving rigid block 


T 

y 


^ Free surface 

Fluid 




Stationary boundary 


Sketch 1. 


impervious boundaries, and they are piaced at sufficiently large 
distances away from the moving boundary to simulate an infinitely 
large region. The moving rigid block is represented by a no-slip 
impervious moving boundary. For the symmetric case, the line of 
symmetry is represented by a free-slip, impervious boundary. At 
the free surface, the boundary condition requires that the tan- 
gential and. normal stresses vanish. 


13 




DIFFERENCE EQUATIONS 


The computational region is divided into rectangular cells 
as shown in figure 2„ The positions of the variables on the rec- 
tangular cells a:t?e shown in sketch 2. The pressure parameter f <j>. 



Sketch 2. 

is evaluated at the cell center* The field variables to be computed 
for each cell are u(i+l/2, j), v(i,j+l/2), and <j>(i,j). The 

cell-centered values of u and v and corner values of uv $re 
evaluated by simple averages. Representative examples are: 

1 

u. . = •=■ (u. ,, . + u. , .) 

1,3 2 1+1/2, 3 z -1/2,3 

(uv) i+l/2, j+1/2 = I {u i+y2,j + U i+1/?, j+l ) 

(V i,j+l/2 + V i+1, j+l/2 ) 

The time derivative quantities are approximated by forward differences 
of the first order 
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/anN D? +: r - d? . 

l,] 

\3t/. • At 

9-tJ 



i+^/2 r. j 


u 


n+1 
i+l/g, j 


At 


n 

^i+1/2 


jlX 


where superscript n+1 denptes the advance tim^ step- All the 
space deriyatives are approximated more accurately by a central 
difference scheme. Representative examples are: 



Ax 



j 



/ 3 2 u \ _ _JL_ 


(u 4 


t u. . . - 2u. . } 


\9xV. . a x 2 i+1 '3 i -1 '! "L'l 

f 3 


/3 2 uv\ 

1 

U^yA . 

^F J 

Ax Ay 


* uv) i+X/2, j+X/2 + (uv) i-l/2, j-X/2 


(uv) i + l/2,j-l/2 - (uv) i-l/2,j + l/2 


] 


In advancing the solution from time step n to n+1, the 
Poisson equation is first solved iteratively by the Helpmann ? s 
methods of successive correction; that is, to sweep along the 
positive x-direction and the positive y-’direction from the origip. 

* ‘ grj 

Ip representing Jig. (12) by finite difference form, the ^ term 
should ideally vanish in order to satisfy the incompressibility 
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condition. However, any iterative scheme will not solve Eq. (12) 
exactly so that the divergence D of each cell will not be zero. 

In order to compensate for this discrepancy, an auxiliary condition 
is imposed in reference (ll) as 

D n+1 = 0 (13) 


which ideally wpuld result in a zero divergence in the advance 
time step n+1. The auxiliary condition generates a self -corrective 
scheme in the computational procedure, so that a relatively course 
convergence criterion can be used for solving the Poisson equatipn 
as tinje advances without the accumulation of error. Imposing the 
condition of Eg, (13) for cell (i,j). 



n+1 

i# 


j 


D. 

1 / 


At 


D n , 
At 


(14) 


the finite difference form of Eq. (12) fqr the pressure parameter 
is 


,n+l 1 

♦i,j = Z 


(<j>, +1 + + $, ■, .*) + — («j». H + 1 + 

^ x 2 i + i/3 1 i/3 .Ay 2 1 /3 + i x /3 i) 


+ R. . 

i/3 


(15) 


where 


Z = 



(16) 
(cont'd. ) 
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R i,j [ (u i + l,j >2 + (u i-l,j> 2 - 2<u i,j )2 ] 


Ay 2 


AxAy 


{v i,j+l )2 + (v i,j-l )2 " 2(V i,j> 


(uV >i+ 1/2, j+1/2 + (uv) i-l/2, j-1/2 


(UV ^i+l/2, j-1/2 " (uv) i-l/2, j+1/2 


(16) 
(concl y d„ } 


D i/1 - 

At 


v 


(D. , , • + D. . . - 2D, . ) 

Ax 2 1+1,3 1-1,3 1/3 


+ -=~ (D. . , _ + D. . . - 2D. . ) 

Ay 2 i/3+l 1,3-1 1/3 


and 


?i,j Ax {u i+l/2,j " u i-l/2,j ) + Ay (v i, j+1/2 ~ V i,j^l/2 ) 


(17) 


After solving the pressure field, the new values of u and v 
for cell (i,j) are computed from the explicit finite difference 
form of Eqs. (2) and (3) evaluated at i+l/2,j and i, j+1/2, 
respectively. The equations for u and v are: 


n+1 

u ‘i+l/2, j 


U i+l/2 


,j + At { E [ (u i,j> 2 - 


■k [ 


(uv) i+l/2, j-1/2 " (i;iv) i+l/2, j+1/2 


(18) 
(cont'd, ) 



, , 1 ,/.n+l ,n+l \ , „ 

g x ‘ Ax vi,j *i+l, j) 


Ax 2 (U i+3/2,j 


* U i-l/2,j “ 2u i+l/2 ) + A „ 2 (u i+l/2 , j+1 


+ U i+l/2, j-1 “ 2u i+l/2 


Ay' 


-] 


(18) 

(concl'd.) 


v? +; !\. ,« = v. . .1 /o + At 

i/ 3+1/2 i,j +1/2 


■ + 


W [ (v i,j > 2 - (v l,j+l )2 ] 

l 

(uv) i-l/ 2 , j+ 1/2 “ < UV) i+l/ 2 , j+l/ 2 ] 

_ , 1 /,n+l , n+1 \ , f 1 , 

g y Ay vi»j ^i/j+l/ V i+lf j+1/2 


(19) 


+ V i-1, j+l/2 - 2v i / j+l/2 ) + Ay2 (V i, j+3/2 


+ v. . , - 2v. . , 1 /0 ) 

if 3~l/2 if J+l/2 


] 


BOUNDARY CONDITIONS 
Moving Boundary 

As the rigid body penetrates into the fluid, the motion of the 
body introduces a nonsteady boundary condition. The position of 
the moving boundary (rigid body) depends on the forces exerted by 
the fluid. An iterative scheme involving the impulse-momentum 
principle and the translation of coordinates is' adopted so that 
at the end of each time step, the Eulerian mesh is translated to 
match the boundaries of the moving rigid body. 


18 



If the pressure and velocity components are known 1 for cell 
(i, j ) adgapent to the boundary, the stress components of the cell 
acting on the boundary by the fluid as shown in sketch 3 are: 
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' r 
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Ay 
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r 



T 



7 T 


'xy 




“w T 


xy 


<P P 


Sketch 3. 


<f>P 


The force components contributed by a single cell are: 


f = CT 'Ay + t Ax - (bpAy 
x x J xy T : 


f y = a y' Ax + . T xy Ay " ^ pAx 


where 
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i _ 


2n 


3u 
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2n 


3v 

3y 
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xy 
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/ 3u 3v \ 

\3y 3x / 


( 20 ) 
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The force components F^ and F experienced by the rigid body, 
treated as constant during one time step, are computed by inte- 
gration of the forces of all the cells adjacent to the moving 
boundary, Equating the change in momentum of the rigid body 
between any time interval n and n+1, to the force acting 
on the body, one observes 

M = Mg + F (21) 


The finite difference form of Eg. (21) for the velocity comppnents 
of the rigid body are: 


u. 


n+f _ 


w 


F n+1 

u n + 

W M 


At 


n+1 


v 


w 


v n + / 
w l 


F n+1 

g + -X 

y 


M 


) 


At 


( 22 ) 


V| t, *1 

g = 0. The values of the force components F and 

X X 

are required as boundary conditions at the beginning of 
each time cycle and must be found by an iterative scheme. The 
displacement of the body at the end of the time cycle, sl^own in 
. sketch 4, is computed from the average velocities as: 


with 

F n+ ^ 

Y 


n 


Ax 


u n+1 + u 
w w 


w 


At 



n+i 

w 


+ V' 


n 

w 


2 


At 


(23) 
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X 



• n+1 n+1 

Since the advanced time velocities u and v as expressed 
by Eq, (22) ijiust also satisfy the momentum equations expressed by 

„ l T 

Eqs. (18) .and (19) , a judicious choice’ of trial values .and 

representing' and respectively are estimated at 

the beginning of the computation. The estimated advance time 

velocities of the body u q ancl v g can be immediately evalu-? 

ated and the reference coordinates are then translated according 

to the, magnitudes expressed by Eq. (23) using u n+ ^ and v n+ \ 

© © 

In solving the finite difference Poisson equation, values of 
<K u, and v outside- the computing region are obtained from the 
appropriate mementum and continuity equation. Although the types . 
of boundary conditions applied will depend on the boundary under 
consideration, the boundary conditions for the case of a boundary 
to the left of the fluid shown in sketch 5 will be presented} » 

The conditions for other boundaries are analogous. 
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Sketch 5. 
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The no-slip and free-slip boundary conditions are the same as 
given in reference 5. 


Free Surface 


The treatment of a free surface cell is adopted according 
to the procedure given by Hirt and Shannon, reference 12. The 
nprmal and tangential stress conditions at the free surface can 
be expressed as 


c£> = tj> a + 2n 


o 3 u . , 

^3u 8v'^ 


3v 

n^ tt — + ni n 
x 3x x y' 

Jy 3 xj 

+ nj 

3y 

- 


„ 8u 
2n m ^ — 
x x Sx 


+ (n m + n m- ) 


/3u 


3v ’ 


(24) 


where the direction cosines of the normal unit vector and the 
corresponding tangential unit vector relative to the surface are 
related as: 



and 


The tangential stress condition is approximately satisfied by 
selecting the velocity components at .the exterior edges of the 
surface cell so that the divergence of the cell vanishes. The 
normal stress condition can be treated accurately provided the 
orientation of the free surface is well defined. At present, the 
free surface defined by marker particles is approximated to lie 
either along the cell boundaries or along the diagonal of a cell. 
The normal unit vector in both of these cases is well defined. 
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COMPUTING METHOD 


The' 'method of this paper calculates a transient problem by 
working through a sequence of small time increments . The results 
of each time cycle are used to define the initial conditions for 
the next cycle. Each cycle itself is divided into the following 
phases; 

1. using a linear interpolation scheme, the trial values of 

hhe force components and F n+ ^ are estimated. Fig. 3 shows 

xe ye 

that the correct choice will lie along the 45° line. Using the 
previous trial and computed values of the forces in the m and 
m - 1 recycles, the trial values for the next recycle will be 
chosen at the point of intersection. For the x-direction 


F n+1 

xe 


p m-l F m _ m-^n-1 

x xe x xe 


F _ p ■ _ f + F 
xe x xe x 


(25) 


where F is the estimated and F is the computed force in the 
xs * x 

m tia recycle. The equation for F^*^ is analogous. 


2 . 1 The estimated velocities of the moving boundary u' 


n+1 


and V 


n+1 


computed from Eq. (22) are used in Eq. (23) for deter- 


mining coordinate translations Ax and Ay . 

w • w 

3.' The pressure ^ ^ for each cell is calculated fpr the 
new coordinates. Using the values of neighboring cells as shown 
in Fig, 4, 


<J>“ . = (1 - Ax) (1 - Ay ) <J>. . + Ax (1 . 

y i,3 w' j w r,3 w V 2 J. 1+1,3 

Ax 

+ % 0- - "T 2 ) h.j+i 


4. Eq. (15) is solved iteratively for the pressure <P'. 

i f 3 +J - 

by successive correction until the convergence criterion is 
satisfied as given in reference 8 as; 
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( 26 ) 



<b' S - 4> s+1 ' 1 

U‘ 
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1/1 1/1 

g h 

* Mj 


th • 

where s means the s iteration and A is the convergence 
criterion which is a predetermined small positive number . 

5. The new velocities of each cell are computed by Eqs. (18) 
and (19) . 

6. The forces F^ +1 and F^ +1 are computed by integrating 
the forces of all the cells adjacent to the moving boundary. 

7. Steps (1) through (6) are recycled or, repeated until the 
force tolerance is satisfied 


,n+l 

. F n+1 

X 

xe 

,n+l 


y 

ye 


< X, 


< A, 


(27) 


where A^ is the force tolerance. 

8. : The apparent viscosity, r] com P u ' te ^ f° r eac ^ cell 

using the cell centered velocity and the stress- strainer ate curve, 
(fig. 1) . 

9, The marker particles are moved to their new positions 
according to the same procedure given in reference 8 . At the 
start of the calculation, the fluid configuration is represented 

by a uniform distribution of four particles per cell. The particles 
define the new free surface. 

10. The cells are reflagged according to the new fluid 
configuration. The next time cycle can immediately begin. 

A listing of the computer program is given in the Appendix. 


STABILITY AND ACCURACY 

All finite difference schemes expressed in explicit fprm for 
initial lvalue problems are governed by some form of stability 


25 



requirements. Roache (ref. 7) indicated that the existing mathe- 
matical theory for numerical solutions of nonlinear partial 
differential equations is still inadequate and that there are no 
rigorous stability analyses-, error estimates, or convergence 
proofs. He wrote (ref. 7) "In computational fluid dynamics, it 
is still necessary to rely heavily on rigorous mathematical 
analysis of simpler, linearized, more or less related problems, 
agd on heuristic reasoning, physical intuition, wind tunnel 
experience, and trial-and-error procedures . " 

The stability conditions recommended by MAC (ref. 8) are: 


C At < 


2Ax Ay 
Ax H* Ay 


(28a) 


and 


2v 


At < 


Ax 2 Ay 2 
Ax 2 + Ay 2 


(28b) 


Hirt (ref. 13) in using a technique applicable to nonlinear, 
equations with variable coefficients, further reqommended the 
following "rule-of-thumb" for the MAC method: 

v < ! At u 2 (28c) 


and 

v < I AxZ (I) < 28d > 

where 

u- is fhe average maximum fluid speed 

3u is the average maximum velocity gradient in the direction 
Dx of the flow, 

It Is necessary in any calculation to consider which one of these 
conditions is the more restrictive, and which will govern the size 
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of At, the time increment per cycle. For calculations involving 
large values of viscosity, the diffusional stability condition. 

Eg. ‘{28b) will bo used as a criteria for choosing the time step. 

At. 

The degree of accuracy desired for a given problem also 
dictates the size of the cells and time steps . In view of the 
iterative, scheme adopted for the moving boundary- fluid interaction 
problem, experience indicated that the time step must be small 
enough to restrict the displacement of the rigid body Ax w and 
Ay^ to one quarter the cell size. The degree of accuracy in 
choosing the convergence criterion X used for the pressure 
iteration has a direct bearing on satisfying the incompressibility 

n t 1 

condition D: ^ = 0. T-he final result will be the force components 
and Fy acting on the rigid body which is the prime feature 
of interest. Generally, if ,Eq. (15) is to be iterated to a high 
level of accuracy by imposing a very small value for X, the 
computational time will be increased many fold for each time step, 
and also as X becomes smaller, the force components approach an 
asymtotic value. Figure 5 shows the relationship, for two cases, 
of X and the vertical accelerative force F^ (expressed in g 
units) for a vertical penetration problem at the end of the first 
■time cycle. Taking g max to be the asymtotic value, the gain in 
accuracy will be less than two percent for X < 2 x 10 -4 for 
case (a) and for X < 2 x lo ~ 5 ■ for case (b) with force tolerance 
X^ = 1 x 10 ~ 4 . Using the computational method in an engineering 
application, it is necessary to establish the relationship between 
the convergence criterion and the feature of interest for each 
case to determine the magnitude of the convergence criterion for 
the computational experiment. 

EXAMPLES OF THE CALCULATIONS AND COMPARISON WITH EXPERIMENTS 

Example I. Low velocity penetration of rigid body into still fluid 

In an attempt to provide experimental corroboration of the 
computational techniques, the impact of a rigid flat block on a 
Still water surface was simulated as an example (ref. 14) . 
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The apparatus shown in figure 6 was developed to experimentally 
observe and record the acceleration, pressure, and displacement of 
a block dropping in water. The water tank itself, a former wipd 
tunnel test section, measures 1.23 m (4 ft) by 1.52 m (5 ft) by 
1.23 m (4 ft) high, and can be sealed to provide for impact tests 
under reduced atmosphere. Viewing ports are provided (as shown) 
for lighting, camera access, and visual monitoring* For these 
testsj water was added to a depth of approximately 48 cm (19 in.) 
and colored with sea marker dye (sodium ffouriscene) to provide 
better visual contrast. 

The' drop model shown in figure 7 was a rectangular structure 
fabricated from ’..95 cm (.375 in.) thick sheet plastic, and measured 
7.6 cm (3 in.) wide by 25.4 cm (10 in.) high, and 61 cm (24 in.) 
in length. The length- to-width ratio was purposely made very high 
tp better approximate the two-dimensional case used in the com*- 
puter study. The minimum dropping weight of the model was 6.24 kg 
(13.75 lb), and provision was made for adding ballast in the form 
of lead shot so that the effects of added mass might be observed. 
The model was held at the desired height above the water surface 
by an electromagnet (shown in figure 7) which could be preposi- 
tioned and released from outside the tank, again tp allow for 
future tests in reduced atmosphere. During a drop, the mpdel was 
guided by teflon-coated steel rods to assure a flat impact on the 
water surface. Coil springs were installed on the bottom of the 
tank to absorb the shock loading resulting from extreme perpetra- 
tion at high .dropping rates. 

f 

Model displacements were recorded by a high-speed camera 
(seen in the foreground of figure 6) operating at 500 frames/sec. 
High-speed color film was used with the necessary lighting 
provided by one quprtz lamp inside the tank and two lamps on the 
outside aimed through viewing ports. 

A reference probe (shown at the left of figure 7) could be 
adjusted vertically to pinpoint the still water surface, and 
vertical displacement of the model could be measured from the film 
by comparing the reference marks on the probe and on the face of., 
the model. A time-code generator triggered an integral timing 
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light in the camera so that time histories of displacement could 
be obtained using a film analyzer. 

A strain-gage pressure transducer of 0 - 103 kPag (0-15 
psig) range was installed inside the model with the sensing face 
flush with the bottom of the model. A ±12 g accelerometer was 
also installed inside the center bottom of the model oriented to 
sense vertical accelerations only, and outputs from these instru- 
ments were taken off the model and out through the top of the tank 
with trailing cables. Both outputs were amplified to increase 
sensitivity at low drop heights and were recorded on a direct- 
write oscillograph. The same time-code generator which triggered 
the camera timing light also, provided a time base for the recorder . 

Tests were conducted over a range of drop heights for two 
model weights. In every case, drop height was established with 
reference to the quiescent water level, and the drop was not 
initiated until the water surface was completely still. Camera 
orientation was checked periodically to insure parallel alinement 
with the surface. The camera and recorder were started about two 
seconds after impact. 

Numerical Calculation 

The arrangement, of 'the computational domain is given in 
sketch 6 . 
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Only one half of the field is necessary for a vertically symmetric 
case. The initial vertical impact velocity and the mass of the 
block were chosen from the experimental values . To test tl*e 
computational method, the mass of the block was chosen small enough 
to allow the hydrostatic fluid pressure to push the block back 
towards the fluid surface. 

Using a square grid, each cell had dimensions 
Ax = Ay = .9524 cm. The fluid depth and width were chosen to 
simulate an infinite medium. The fluid depth measured from the 
face of the block was represented by 24 cells (23 cm) , and the 
width of the fluid for a half field was represented by 60 cells 
(57 cm) . The height of the block was chosen to assure that at full 
penetration, the block would not be submerged in the fiyiid. Dimenu 
sions for a half-block were 4 cells (3.8 cm) wide x 34 cells 
(32.4 cm) long. The kinematic viscosity of the fluid was 
0.1 cm 2 /sec. The mass used for the two-dimensional case was 
calculated from the model weight per centimeter of length. 

Comparison of Results 

A comparison of computational and experimental results is 
summarized in table 1 and presented in figures 8(a), (b) , and (c; 
to illustrate the effects of changes in drop height and in body 
massi In figure 8(a), displacement, acceleration, and pressure 
time histories are shown for an experimental drop height of 1 cm 
using a model weight of 6.24 kg (13.75 lb). Analysis of the film 
used to measure displacement indicated the initial impact velocity, 
v , to be 33 cm/sec within an accuracy of ±5 percent. 7?he 
computations were conducted for v q = 30 cm/sec, with convergence 
criterion, X = 2 x 10 and force tolerance = 1 x 10'\ A 
time increment of At = .003 sec was used for the first five 
cycles, and was reduced to .002 sec thereafter as the velocity of 
the model increased. Computations were terminated at t = .411 sec 
after the model had attained maximum penetration. As can be seen 
in figure 8 (a) , the agreement between experiment and computation 
is reasonably good for the water penetration phese. Calculations 
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were also performed with the space grids, reduced by one-half with 
Ax = Ay = .476 cm and essentially the same results were obtained up 
to time t = .22 seconds after which the computational method 
became unstable. 

The effects of increasing drop height to 2 cm are shown in 
figure 8 (b) where much the same trends are noted as for the 1 cm 
case. In this instance, initial impact velocity for the experiment 
was found to be 47 cm/sec ± 5 percent, so for the computation 
v Q = 50 cm/sec was used. The convergence criterion, force 
tolerance, and ‘time increments were the same as for figure 8(a), 
and computation was terminated at t - .395 sec. 

In figure 8 (c) , the mass of the model was increased to 
153.5 cm,, and the drop height was 1 cm. At this weight, the 
experimental initial impact velocity was found to be 30 cm/sec,, 
and v q = 30 cm/sec was also used in the computations, i^gain,' 
the time increment,’ At, was .003 sec for the first five cycles, 
and was reduced to .002 sec thereafter. At time t = .155 sec, 
the displacement of the block exceeded one-quarter cell size, so. 
the time increment was reduced 25 percent to .0015 sec, and compu- 
tations were terminated at t - .21 AS sec. Again, reasonable 
agreement is noted between experiment and computation. However, 
it is felt that the agreement is sufficiently good to validate 
the program and give confidence to the computational techniques 
employed . 

As previously suggested, the computational method not only 
gives a solution for force and- displacement of the block during 
impact and entry, but computes displacement, velocity, and pressure 
throughout the flow field. The computed behavior of the flow 
field cannot readily be examined experimentally, but the agreement 
between computed and experimental block behavior lends credence 
to the predicted fluid dynamics. As an example of how the compu- 
tations may be used to examine in detail the evolution of the 
fluid dyiiamics during water entry, figures 9(a) and (b) presents 
particle' and velocity plots for an initial block impact velocity r 
v q - 50 cm/sec. 
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The fluid pressure developed under the block may also be of 
interest in certain applications. Figure 10 shows how the excess 
hydrostatic pressure developed along the centerline of block 
varies with fluid depth as measured from the bottom of the block 
as it contacts the free surface. Two initial impact velocities 
are shown, v q = 30 cm/sec and v q = 50 cm/sec, and the figure 
indicates that the size of the fluid field originally chosen was 
adequate to simulate an infinite medium, since the excess hydro- 
static pressure is insignificantly small. 

Example II, High velocity impact of a. rigid block into a 
very viscous non-Newtonian fluid. 

The arrangement of the computational domain is given in 
sketch 1. This numerical example simulated the case of an air- 
craft wheel landing on soil runway with a locked-wheel braking 
condition and zero lift. The dimensipns and mass of the rigid 
block represented u two-dimensional equivalent of the case given 
in the figure 21 of reference, 2. The stress and strain-rate curve 
for the viscous non-Newtonian fluid, simulated a soft clay soil 
with water-content of 33 percent (fig. lid of ref. 6) as shown in 
figure 1. 

Using a square grid, each cell had dimensions Ax = Ay = 2 cm. 
The fluid depth measured from the bottom face of the block was 
represented by 10 cells (20 cm) and the width of the fluid was 
represented by 29 cells (58 cm) . The height of the block was 
represented by 4 cells (8 cm) and the width represented by 7, 9 
and 11 cells (14, 18 and 22 cm) to simulate various contact length. 
The mass per centimeter of length was 90,700 gm. 

At = 1 x 10 “ 5 sec was chosen using the diffusional stability con- 
dition and the maximum flow parameter, t] , obtained from the 

a PP 

initial secant slope of the stress and strain-rate curve. This 
small time step will impose a severe restriction on the practical 
gse of the method because of the large amount of computer time ■ 
required for a given problem in which the time history of the force 
and displacements of the rigid block would be long. Only sb ort 
calculations were made for this example, for t. = 0 to 
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t = 1 x 10 4 sec, using the same convergence qriterion and force 
tolerance of example I. 

Figure 11 shows the time history for the horizontal and 
vertical forces experienced by the block for various contact lengths 
In all cases the initial horizontal impact velocity, u = 500 cm 
per second, vertical velocity v q = 0 and the gravitational force 
are instantly applied. The horizontal force, developed by the no- 
slip boundary condition, is proportional to the contact length- 
The magnitude of the horizontal force decreased gradually with time, 
whereas the magnitude of the vertical force oscillated upwards. 

The horizontal and vertical displacements would be very small. 

Figure 12 indicates the displacements at the end of 
t — 1 x 10” 4 sec. As expected, the horizontal and vertical dis- 
placements decreased as the contact length increased. 

The yariation of forces with initial horizontal impact 
velocity was calculated for case with contact length = 14 cm. 

Figure 13 shows the effects of increasing initial horizontal, impact 
velocity for the time interval between t = 3 x lo ~ 5 sec and 
t = 1 x 10 sec. As expected, the vertical force remained con- 
stant as the horizontal velocity increased. The horizontal force 
continued to increase as the horizontal velocity increased. Since 
the sinkage of the block (vertical displacement) was extremely 
small, the horizontal force experienced by the block was contri- 
buted mainly by the viscous tractive force and resistance asso- 
ciated with the inertia effect or the "bow-wave" did not materialize 
The viscous tractive force increases linearly with velocity. 

CONCLUDING REMARKS 

The viability of the computational method for studying the 
low velocity penetration of a flat rigid block on a fluid medium 
has bqen demonstrated to agree well with experimentation. The 
method is capable of providing the complete pressure and velocity 
fields, and the visual display of the history of fluid 
qonf iguration . 
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Some limitations inherent in the method must be mentioned. 
First, only non-turbulent flows are considered in the model. At 
high speeds of impact, the turbulent effects may not be ignored 
and some implementation for turbulent simulation may have ho be 
incorporated in the method. Second, if the main feature of inter- 
est is centered at the short time impact phase, the choice of 
time increment may require some computational experimentation. 

This limitation will not be critical if the main feature of inter- 
est is centered at the maximum penetration of the impact body. 

In the comparative test cases by the numerical calculation and 
by experiment, the viscous effect is relatively insignificant 
since the viscosity of water is small, and the forces acting on 
the body are contributed predominantly by the pressure. 

The computational method has not been tested for the case of. 
fluid-solid interaction involving a highly yiscous non-Newt:onian 
medium. The apparatus to provide experimental corroboration of 
the computational technique for the case of example II would be 
complex and no experimental study was made for this case. 
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sr x=total forces i k x-niitFcmiN 

c_ 

5FXG=“sn»Arr of r j n al hiirizdniai whffl forces ...... 

r. 

Y-roui i in y-i>i iif* 

c 

$F y ti— > XT I FA r ( l)F 1 INAL VERTICAL VHFFI rnPCFS 

r. 

rignaxi i.ji= nornai sri'ess in /-direction 

r <■. i r . .iaym . 1 1 = MiPMAi athfss in Y-i'inrCTinN 

£ 

si.,mah i . j » -M4 x i mijm pRicrpt.r stress of cacf cm 

c 

ST-HAYsCayHIUM STRAIN FNC.CUMTLPFn nil" INC. 0 YCt F 

C 

ST/\^ (l./SkC) DATA PUT^TS 

c 

SI !*■ S= S T3£ SS I DYA.FS/ I CMP*/ 1 1 OATS POINTS 

c 

SUM = SllP''AT10N Of ABSOLUTE VALUi CE RII.J) 

c 

t - = 

c 

TAP FI A )- OUTPUT STORAGE OF VFLCE.tTI^S 

Q 

TAPEI91- lur^ur STORAGE OF STRFSSFS 

£_ 

TAPE1I0I- OUTPUT STORAGE up PAPIICtF C OOPO ( 6 A TES 

c 

TAP-Ill)- T c P D nPAPY storage of velocities of PREVIOUS TINE STEP 

r 

TAPE! 131- OUTPUT STOP ACT OF PMI.FTA.f.FF 

c 

T10EU4)- TERPOPARY STCPAGE OF PA«T 1CLE COORDINATES 

c 

TAPE! 15)- TFPPtlPAPY STORAGE OF PARTICLE CCGR01NATFS 

c 

TAPEf I d) - TEMPORARY STORAGE OF VC10CITIES FOR RESTARTING 


TARTY! 1, 1 1= SHE \R STPrss ON X-Y PI ANF 

r. 

T0P= TI«F FOR CFI.L PRINT 

c 

T">-;=TI4F 1 A TA ^rnkACP 

c 

TH= TAI I . Jl’ANGLE OETWEEN POSITIVE X-OIRECTION AND SIGHAHI.J) 

c 

TI =TIHF TO Oil I T 

C 

TPP= T r . 1= -HR PAPTICLP PRT^f^ 

r. 

TS^= T I *‘F CPR STRESS PPTfJT 

c 

ill ! » J ) = X - V - L lie I T Y OF CFLL 

c 

UG==STI^AT= [)F FINAL HORIZONTAL irJHEEL VELOCITY 

c 

.niii.jiM-'/cinriTY at time stcp m+l 

r. 

tJft = X-VELOC ITY OF .'HEFL AT TIME' STEP N 

c 

VI I . l)=Y-V=LOr. ITY OF TELL 

r 

VGs5STIMAT= '"‘F Ff^M VERTICAL lit'rFl V^LOCTTY 

c 

VNI 1 1 > J ) = v— V** LOC I TY AT T l Me STFP N+l 

c 

V .'a Y-VFf or I TV OF PME EL AT TIME STFP N 

c 

F'M=.MASS OF JHEEL PER UNIT l.IDTH 

c 

VM=SOUARF IF DISTANCE PROM ANY POINT TO WHEEL CENTER 


UXsI-r.nOPf IN4TE OF UHEEL CEMER AT TI"E STEP N 

C 

•SY=Y-r.n 1RI INATF TF WMF FI CFf TFR AT TIME STFP A 

r 

xn )=x^rnr->fiT*JA ft np Ctu rF^r^F ... . 

C 

XK f K) = X-C'-1P0INATF OF PARTICLE K 

c 

XO^ INITIAL X-COOROINATE OF FIRST PARTICLF 

c 

X PL 1 1 ) =X-C JORO I NATE OF CELL RIGHT BOUNDARY 

c 

XB*X-CnORCl')ATE Of RIGHT BOUNDARY 

c 

<S r PAR = CO , nre , l FC° NO. OF GETTING COLUMN OF PARTICLES 

c 

YIJ) s Y-C0f- , O INATF OF CELL CEMTFP 

c 

YK(K)=Y-CC GROIN ATE OF PAPTICLF K 

c 

YO=t«lTIAL Y-COOROINATE OF FIRST PARTICLE 

r. 

Y<>L< J)=Y-CjOPD!MArE OF CELL LOWE” BOUMDARY 

c 

YSF°AR=COL‘ ITER FOR NO. CF GETTIM ROI! OF PARTICLES 

c 

YT=Y~CnC^r,TMATE OP IOWER BOlJNDAP Y 

c 

OI Mr\S ION UI2 1.171 ,V(3 l, 17 1, PHI I 31 . 171 .»! 31, 17),F (31,17) , 

*FF1 3 I. 17 1.Af.( 31. L7).FTA( 31. 17} .FK(31 ,17) , FCI3 1 ,.L7_> , 

AS1G* AY (31 ,17) ,SIG‘AAY<3i. 17 ) , T WXY t 3 L . 17) , X<( 1020 ) , YK 1 IC20) , 

fX( 31). YM 7>.yiL( 31 ) .YPL1 17),NA.''b< 10) .STRANt 56) ,ST3FS<56, 1 1, 

*STRAI\'t 1) . STP=SS( 1, l ),0FLY156,1 ),H15S> ,S2< 56, 1 ) , $3 1 56 , 1 > , 

=*TXM 10 30), TYk ( 10201 .!■( 31 . 17) .DATA!™) .TUI 1020). TV! 1020 1»OAT3(301 

EJUlVALP^Cc l «*TU> ,*(YK,TV > 

INTtGFR F,FE.PC,F$ 

C _ 

c 

REGION 11 GENERAL SETUP 

C 

4000 COnTI « dF 

R £ AO (5.11 ('!.'•£( 1 1 . 1=1,8) 

IFIE'F,5)99C9,4031 

4001 COM I'M- 

RE Al (4,2) It>AC, JBAP.NIS.NIF.NJS.NJF 

READ (6,3) £TV>, FKO,GX,GY,XR,YT,RW.UX.KY. WM , UG , Vh , NX, NY, XO, YO ,DXK , 

*l)YK.r*.rlY,< 1. iTCP.f T P°. 0 TSD , ALPHAIJ . ALPH AV , CPNV,G APA F, SK'S, [WHEEL, 

*1FSY“,FS. [HA<,.VI 

RFAFt a .11) -i-IPTS.Mt'CVS, UfAX,NS,NCVS,MS, IW, ITG, EPS 

prA0(N.?0) I'rPAYM I l.STRESM ) .1=1, NS1 

R6AO(5,jo) I,TL,NCS,NSP»NPP 

l FnR^Ar(*i.M) 

000«> 

3 FCR‘Mr (l.= 12.4/6 = 12. 4/2 112. 4E12.4/6C 12. 4/5E 12. 4,4 1 4.F4.C) 

0004 FOR • Aid'll// 1 <.8\a.//.6X ,2Ht =, 1 l2.frX,2HJ=. t 12, 3X , 5HETA0=, E 1 2. 4, 4X. 

*4 <f ^*"'= » ‘ l ?*4t y \ t **IHX= t r l2*4*5Xt3l’GY = » ,: l 4 * / T 5Xt i l ?Xrt=if 1?.4» 5 X » 

*3KYr = ,';i2.4,5<.3'|R«=.' : 12.R,5X,iM5X= t !: 12.4,5X,3HJY=,E12.4,‘)X. JHWM=, 

*112. 4, / ,5Y, IF IU-, c l?.4,5X,3HVW=,r 12.4, 5* . 3HMX=,I 1 2 , 5X , 3HRY=_, 112, 

•> 5X , «,■>' l=,r 12 , 5 S, ,HYl' = , r 1 / . 4 , / , 4 X , 3 XK -- , t 12.4, AY,4F(YK=,F!2.h, 

Y5X, M- 'X = ,F 12. ,, iX. 3HflY = . F 12-4 ,SX, M UT= ,f l?« 4, 3 X , 5H0TCP=, F 1 2 . 4 ,/ , 

*3< 1 4il<'TOP = ,si2. 3x .5Hl'TSP=. 312.4 , 1 X . 7HALPHAU* . L I 2.4, IX, 7HALPHAV=, 
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*j:.l2.3i-iSiJ5 , *C;_'IVSjFl?,^,^i6HCAMAr = l Pt2.4,/,3X I 5HHHClS=fri?. 

* j)iT=.i- iz.i.m. jurt-.r !2.4,'.x,4iiNrs=,i i - a. x nsp= » j i ? , _______ 

*4x, >».')»!>-, I I IX, /HlWMti t»» 1 12 , ?X, <,ij_[l. S.VX? ,J,12 tJXj 2 PAS^jJ. 1 2. i — _ 

12 X , « 1 IJ:f '_X = , li 2 r 5 X , 3 1 ! K L= , U 2 . 4 , / /V ) 

S I I'..’ ■'ATI /////■■•»■ I <|I Till ST' P T-.t l?.4.2y. UH TOTAI M1WHFK 'if : 

«1M « 1 1 OLSZj I 1 ;.J 

00 06 I -UpAT HAS flCWF rHF<lJ._TW2 sr AGE s , 2X , 

* J r L .-t’-.a.LA ; F f £ T F I'll-S, /-x.lKii MI =., lii?.X x.Z MJ r ,.12,_2X,.2JiF^L2j 

S2X jj ,_Li, 2Xi3;{SU:j.U2, JJ 

0.007 F.'-UElVt 1 1 = ,J Z.^A^-ZlU^qz, ZZ*.l 11S.LU.JJjl,.£JL 2 -,2X *ZhJU J_aJ_l= ^ 

11- : L2.4) 

omo riN"tn///,;',;3'OM i, Ji r.H gr k at t ms t=, e l 2 *4., /j 

g.0. yj I UZ&* i:!C7.C.LE_J^_,iUZ i .i J 2X,3l!ia£GAf-LA.L_ L f:j-a..2 J 21 

0010 hlMiXAJJ vx . • I T - . F 1 7.4.<i X, 7H 1C = - I 12.4X, 7h a.AV3,F. 1 2.4, 4X. 

tiii SJX±*ZiC-LL,2jJ , 

!LCUl_£l!idiiAT 1 -x.3< r=. El 2.4 .4X. 7* ^L^iil2,At-4 X.,_7 h S.fy=, t jiz , <U 

001? rfl- J -^ATt//.0V.. > -jO r CMPLI TAT H 1 ^ FMPei> AT TfJF Ts.rio,?./! 

0013 FOPdATt 4X.:-I r-g.Fl ?. W.X,7Ht -' B;10HU=,l 12.B,-'tX, 7 m.nRnFV = «H2.3, 

* 4X.7T‘ tl'U' .11? ) ~~ 

0014 F t |p- -« AT ( > X. - WE-CYCLE TC IWP PO VF VEtOClTY E ST IMATE AT T*,F12.4, 

*/ , 2Y,40'Jw'l = .E16.0.5X,4KVWM =.~i4.a i_6 X . 3HIJ |0* ,_E’l6 . 8_, C.X_,j H VG = .E16.8, 

*/ , ? <« 4H5FV=.-; 10. 9 , OX ,-tHSF Y^.F Ife - 0,?X^7Fi r RR0PU = . 0 16.3, 

*7UF"7m iV= .Flr.ri. >X. THIO. 1 10) I _ . 


00 13 FQJ ''AT (/■?<■ = « >HH--H»I 70NTA 1 HHFFl V E KNOTTY 7F*7n AT T-.. F 1 2 . 4. /, 

*2X . 1 5 u C A E f.lurF fl S FXs.Fl?.41 I 

0015 FF|Q J ''ATl‘.X.?3-\0.IJSTr-i- TO'E V Aftl API FS. //-.4X . 3F0T= , FI 2.4, ?X. SHOTCP 






















onn ono ho 


























ono nbn nob 


K<JAft_;NX*NY ; J 

jp /JtOPO. ; 

'1 P R b- K M ,\ I< - J PA R± ljl -L0 

MpAR|»«ipftatI ~~ t 

"V r* — ‘ j in” ' 








Ut 1 j J)_=1JW 

y ( i . j i =Vrt . 

LL( hsv-i.eQ.i i on rn na 

LL i i^.Jl =5 

»1 1H-1 . J1-1U 

iiLLdiJlj' uh 

VJL * v a 

clip co m f iuu- 

rp (~n-~R»7r*ox*(x( i i-oxi-wxt*? i .if . mw <-»?n r.o rn 121 

IF t ru+? . »nx» 1X( 1 )-UX 1 H)X»*?-.>.»»YaiVIO )-«Y >4CV»«7).I F.(BO»eZl 1 

*r.n T'l 1?? 

FlIh i .ji-; 

tr 1 IFS Y -i. eo. 1 1 no re 121 

FFl 1 1-1, J 1 - 2 

00 rn 1a 

0122 fc t r » 1 .J)=2 . . 

F£i m . 1-1 ) = 1 . 

n i irsv^.go. 1 1 r,n to 121 

fc ( t“ -i . Ji =2 

FE t'lH-l . J-l 1 

121 CUNTINOS 

1 1 1= 

126 crMir-ij:: 

HM=(FXFRO*COS( i. 1416/4. 0H/DX+3.0 0 

N=N0-1 

NV= ( WY -FRO 1 /DY *2 .0 

m=i <-l 

00 121 l-rl,Q«-l 

1)0 174 J = .V.F.M 

lifts 1X1 I)-0X)*y2+IY<J)-WY)*s2 

[FIWP.GT. (Rrt«v2) 1 00 Tfl 124 

F1I.JH5 

Ul T . J 1 = 00 

vi i . j) = vu' 

1F( 1 FSY'1 . EQ . 1 I GO TO 139 

FI IB. J) =5 

01 I M- 1 . J ) =0W 

01 1B,J)=UH 

VI IB. J 1 = Vd 

0130 CONTINUE 

IF 1 ( f2.Q--0Y»l Y( ■J>-V'V)»0Y*’C21.1F. l 0 «<-42) ) GO TO 124 

1 F t 1 2.0<-QX*( XI 1 l-ViXU-OX 4*2+2. 0*DY?1 Y1 J ) ->.'Y ) +CV4*2 1 .IE . t P0« *2 1 ) 

*GQ TO 125 _ ___ 


FFl 1.J4-1 ) = 2 

IF! IFSY-I.EO.l 1 GO TO 124 

FE ( H'.J + 1)=? 

GO TO 124 

0125 FFl 1.0 + 11=2 

FFl 1-1. J+l 1=1 

IFllFSVi.EQ.ll GO TO 124 

FFl IF. J+l 1=2 

Ftl 1“+1,J+1)=1 

124 CCNT1.MJE 

Irt= IB-1 

' 127 CON T PIPE 

0113 COKTIM'JF 


C SET-0 P CB ADJACENT TO 8ND CELLS 


IFINCS.Nc.Ol GO TO J30 

DO 14 1 1=2,1FAF1 

IF 1 FII.2I .EQ. SI GO TO 143 

FMI ,21 = 2 

0143 lFinirJ i lAKll.FQ.51 GOTO 1-1 

FE 1 1,004211=2 

0141 CPMTHlOE 

III' 14.2 0=2. JE AR1 

lF(FI2i J1.EQ.51 GO TO 144 

FE 1 2,01=2 

0144 IF IPltW 1, J 1 -FQ.S) GO TO 142 

FE 1 1 FAR I . J 1 =2 

014? rPNTI.VQF 

0 130 CCNTIf.-OE 


C Sl.T UP SUR CEIL AT INITIAL CQttDl T10NS 


1F1NCS.MF.0) OH TO 163 
J=FS 

on l?i 1=2.111/011 

IFIFFH, Jl.F'J.2) GO TO 131 

Fi l.J 1-3 

111 CUNT HU JF 

lFllfsV1.FiJ.il r,n TU 140 

r i < , ) ) ~* i 

111 2 » J ) = 0.0 

V12 .01 -o.o 
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01 40 C fHIIfr’>lF 

FJ I l'V> 1 . jT =T 

U 1 1 1* A-' t , J )_=0 ._0 

y.i.i i>a.( 

01 a a n‘Mfi.M'ir 


C rSTAOLtS lI INITIAL H YDKPST A T 1C P PtSSHHF 

C 

irtNC s .oF.o) on t7i ti4 

N^FStl 

oo i n-i j~M , j-* as l 

0 I 10 0 1 -?, 1 T AX 1 

PHI ( I . J r.CY-IJI-IWY + RW+.Sum I«GY 

0109 CCNIlU'Jt 

0114 CUNTIHOF i 


C READ nt IJ) .</(!J) FROM TAPFI 8 ) AN P STUPE THESE 09 MAIN STORAGE 


if mcs . eo.oi on rn 15 / 

0104 SFADtf.lDATA ,U,V 

iFtr.'iF,ai 104, io8 

oios coktim-j= 

nis=n»TAi n 

Kii I ? ) 

SI>XSF = PATA(3) 

SOYSFeJAT A (4J 

XSr PAR-OATA! 5) 

YSF P AS =04 T A ( 6 ) 

SFX=9ATA( 7) 

SPY = HATA ( 3> 

Vri^DATAl 9 ) 

U4=CATA t 10 I 

ALPFAU=DATA( 1 1 ) 

ALPHAVsQATflt 1 ? ) 

nOTfJ=OATAl 1 31 

QOTSP-nATAU*.) 

OPT PP=0 A T A 1 15) 

OX SF=PAT A 1 ?5 1 

DYSF^OATA 1 ?4) 

ERKQRU=PATAt ?7) 

EHPQPV=r>ATA( ?A ) 

IMDATA129) 

SU‘'3=r>ATAt30) 


SFXG=RATA( 31 1 

SFYG=OAT A [ 32 1 

IF 1 A0S( TJS-T1 ,LT. ( DOrCP/2. 1 1 GC TO 904 

W 3 I T c 16t371 

GU TO 901 

0904 PEAC1S) 

IF( EOF, 31 155,904 


C REAP FIKAL SET DATA OF PH 1 tl J ) ,F ( I J 1 , FE 1 t J ) , AND ETA t 1 J 1 FROM 

C TAPEH31 AND STPRF THESE CN MAIN STPRAGF 


0155 SEA0113) O ATA,PHI , £TA,F,FE 

1FIE0F, 131153, 905 ; 

0905 T05— OAT All) 

01)TCP=0ATA1 131 - 

1 F ( ABStTOS-T 1.LT.1C0TCP/2.11 GO TO 903 

WSIT’;ifa,57) 

GO TO 901 

0903 REAM 131 

1 F t EOF ,13) 156,903 


C ESTABLISH COUNTERS FOR PARTICLE COORD. CALCULATIONS 


0156 SEAC114I0ATA 

1 F ( FI'F. 14 ) ISA, 153 

01 SB CONTINUE 

Ti:S=DATA( i 1 

KB AS -PAT '3(21 

OCTPP = ')Ant 141 

1 F 1 AOSI TPS-T I .LT. (ODTPP/2. 1 1 GO TP 902 
HP 1 TF< 6,571 

on rn 901 

090? MPAR-K3AP / I0?0 

HP9p-KeAft-MPrtR*1020 

MPAS 1-9PAH-H 


C PRINT INPUT DATA EPS VEP1FICATI0N 


0157 HUTF ( •, ,/,] ( N A M L ( 1 )_, I = 1 ,Bl , nAK.JRAR, FT AO, FKC . GX . GY , XR , YT ,S W ,HX , 
_ ... *riYr.H'l»iH» VV.,MX f NY,XO,Vo,DXK,PYK, OX , OY , PT , PTC P , PTPP , 0 T SP, ALPH AU , 

*A1 Pil W,r'r]7'iv,l'iAHA'-: tSHPS.Y. 1L,NCS,NSP,NPP. 1WHFLL .1 FSYH.FS, 1THAX.WI 

HR I 1 F 1 b ^4 0_3 2J 

WS1 II (A, 40 ?~7 ) HNPT S.HiSC VS t MMA X. NS. NC VS f PS t 1U, 1 TG, EPS 



WrtLTlLJAiJ^flJ __ 

HR 1 in 6,400V )1 I , S TRAM! I > , S H?c SI H . 1 ~ I , NS 1 


I STORE INITIAL VALUES ON TAPES 


LO_NCS .NE . 01 GO n) 1/1 

wm rn & j/al 

A S SJ ' . K 1 70 TO KKf T 

TPP-'l. 

(~,n n< alt 

ni 70 A5S1CM Ifl TO KRFT 

rci’^o. 

r.n to ah 2 

0171 C0N71MUF 


£ (7CG10N 30 RIIJ) CA1CULAT10M 


300 T= T-t-DT 

CM.L SECOrJOlHT IME1 
HRITfcl6»9)T.3riHE 


C CALCULATE L) FOP COMPUTING THE SOURCE TERN 


004014 J=2,J0ARI 

009034 T-Z.IEWUM 

on «ji =q. 

in f< r . j ) A'-io.Ft i . j i .op. 3) c-o rn 4034 

D ( i . .1 ) = ( 1 ill i . .1 > -u< f -1 . j ) i/iix 1 + t ivn . j)-vi 1 . j- 1 ))/ny ) 

4034 CONTIMOr 


C CALCULATION OF THE SOURCE T£OM 


SUHftaQ. 

00 301 . J3AR1 

00 1=?, in fl R! 

3 1 I 1 J ) =0. 

IEIF1UJ) .WE . 2 I GO TQ 301 

U5~= 1 IJ< [. 1 l-HJl [-1 . J ) ) **Z 

vs=i v( 1,1 i+vt 1 , j-i ) 1**2 

uvoRnui 1 , j+i 1 *ut t « j> >ffi v < i-n, j)»y u i J 1 1 

uvtl=<u( t -1 . j mil i-i . j-n m vi 1 , j-i ) +vi 1- 1 , j- 1 > > 

UVT°n 01 I »J 1+ 11 I , J-l ) ) *1 V( 1 t 1. J-l 1+V (T . J-l) ) 

UVBL-HH1-1 t J»1 )-HJ 1 1-1 , J ) >*t VI Ll J 1 + V< M . J ) ) 

DLEFT=Q(I-1.J) 

DRIGHT-DI 1 + 1, J) 

DA80VE=0( 1 1 J-l ) 

D8EL0W=I)( i. J-H 1 

1E(FE( 1 t J 1 . NE . 2 ) GC rn 319 

I ft FI!- 1, J l.E 3.3 .Hk.FI 1-1 . J1 .SC.S.CP.Fd-l, J) . EQ.6) GO TO 30Z 

ULS-1UU-I. JH-UII-2, JH*«Z 

GO' T.1 303 

030? ULS-US 

DL E F T=0'( 1 ,J1 

TFt Ft 1-1, J> .E0.5 ) GO TO 304 

UVPL=0. 

U'/XL ^ . _ : 

GO T n 30? 

0304 UVPL=LMfty t )S4 . 

UVTL»UW*VW*4-. 

ICIFF! 1-1 . J1.NE.1I GO T r. 303 

1F(F( t-l.J-D.nlE.Sl GO TO 305 

ORELru-OflELCUnUt i-i, J+l l-UWl/DX 

GO TO 30 3 

0335 Drt8'.)VF=rUBnVS+(Ul 1-1, J-l 1-UW l/DX 

0303 IKFim.Jl.iO.l.OB.HIil.Jl .E0.5I GO TO 306 

U3S = U1( H-1 , J) »|J|1 , J) [ 

_ GO in 10 7 

0306 u:is=us . 

DR IGHT°OI 1 . J) 

ififi in ■ j) .h ).5i go rn 308 

UVlTR-O. 

IIVT'I=0. 

on TQ 30 7 

0300 DVflR-LM*VH»4. 

uvT><=uo*yw»'«. 

IEIFE1 1*1 .Jl.'lE.t ) GO TO 307 

lHFn + liJ-l).NE.5) r,n TEl 300 

imSl ’14^00 FLGk^UW-lH l,J>l))/f.y 

GO TU 307 

03 30 l)AUOVr=l)AUntfL + (UM-UU.J-in/t>X 

0107 If (K 1, J- H.f ). l.nK.F I 1 ,J-1I.E0.!)1 GU rn 310 

v*s=(v( iij-in-yn, j-2ifo»2 

GO TO 311 

~OltO V4 s — VS 

r uHOvt =<M i ,ji 

iririi,i-ii.h).5) gc rn 3iz 
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IIYTI.M, 

ijyjK_=o . 

wi ,_Tji. ju 

O J.L? '■'ll VTlJ=J.l<j»V.w»4 J _ — _ 

wvT»=nyt*V'i»<t. 

11-IK.I I ..I-H ."IF. L 1 co rn 3 1 1 

ipim i — i.~j- 1 1 . nj .>) r,n t ii 3i 3 

fjjy.nHT 3 . 0 R iqhl + L/J i± 1 ?J - 1 )-vv< )/nr 

fi'i to j 1 1 

— 0.3 13 umiiOL 1 : rr -j/.wiv.Q. 1 ; 

mi ) mni f IH) 1 . Of’. u. on T n 314 

V RSMV t 1..H-1 )*■'/< I.J) )**2 

till 111 113 

0314 V0S=VS 

OltCI nn=tj ( t, J 1 

__1 J- IFI1. J + l ) .F U.S I on TO 316 

Jvnt=n. 

nvii^o. ~ 

op TO 313 

0316 tlVflL=U/« ; VH» 4. 

UVPR^UW4V^*~4. 

[FIFM I, J»1 I .NF.l I on TC 315 

II I f f 1-KJ+ll .NF.5) GO Til 31B 

1)0 3LGH1 = OR I GOT f ( Vll - V 1 I -t -1 < J I 1 / ny 

r,n Tn 315 

0310 [)LEFT-OLFFT+(VH-V( 1-lt Hl/DY 

Gil TO 315 

0319 T MFMl-l ,J-1 ) .N£. If CO TO 320 

JVTLstNfcVyH-9. 

DABnvF-'MmV£HU( t-L.J-ll-HH t-l. J )-(2.*lrfH/rX 

DL!iFT= JLOFT-t-CVn-l. J-l 1 + V( I , J- 1 ) - ( 2 . * VWl 1 /OY 

GO TO 323 

0320 IFIFM I-l.J-H ) .NF. 1 I CO TO 3Z1 

uveL=n9*vvJF4. 

01 f F T = -)L-FTM I 2.*tfW)-Vf I.J )- VC l-l.J) ) /OY 

D9fcLrw=IJOtLnw+(Ul t-l .J+l )-HM t~l » J )— 1 2 . *UW 1 1 / EX 

GO TO 323 

0321 iriFrm-i.j-n.HE.n go to 322 

UVT)iiiU.JtVW<=» . 

flAPCVE^ABHVE + t (?.»IJHI-U( 1 , j 1 — u { 1 - j — 1 ) )/nx 

00 1C-HT = I)I' IGMT + I VI Hli.l-11-nfl I « J-11-12.*VVi) ) /OY 

Gn TP. 323 

0322 IF1F C (H-1«JMI.NE.1) GO TO 323 

UVBP.-U9*VWft4. 

1'R1GHT=ijR1GHH-I 12 .*V2) -V U , J >-V( 1-H, J) )/DY 

D Q ELOl-OG £L0». +1 ( 2 . <-Ua' 1 -U 1 I . J 1-U 1 i , J* 1 1 ) /l) X 

0323 URSMUl t*l,J>»Ul 1 .J) t-»2 

ULS=(U1 1-l.JH-m 1-2.J1 l<=*2 

V5S=H V< 1 , J»1 1 *VH , Jl )*=*? 

VAS=(VM.J-l>+VU.J-2 ))<■<•? 

0315 HI ,J) = lUBS-HJLS-(2.=fUS>l/(4.*PXS) 

*4-<V6S+VAS-(?.-i/S) )/(4.*DYSI 

*MUV9f'4-UVTL-UVTP.-OVBl.)/<?.4QX*m-Dl I » J 1/OT 

*-ET A 1 1 . J] *1 DO IGHT ■>0L£FT-(2.*0(tt J 1 1 1 /P XS 

<--ETAIl .J>*{08£UH + iOA6PVE-(2.*0<I.J) D/OYS 

SUF R— S'JMF. t AOS t R 1 l.JII 

0301 CONTINUE . 


C REGION 40— A FREE — SURFACE PRESSURE CQPRFCTICN 


DO 451 J=2..IBAR1 

nn 451 1=2*1 0AR1 

IF1FI t,J) ,NE . 3 I GO TO 451 

_ PHI ( 1, JlaQ. 

IF(F(1.,J-1>.N'F.4I GC TCI 45 2 

IF( Efl-1 , J) ,NF .4) PC TP 53 3 

IKFf I + I.J).E J.4.-' , R.H I, J»l) «E0.4) GO TO 951 

phi < i, j)=p. 3 *fta( t , ji*( me t, j+n+m t-i.j+n-ut i.j>-u< 1-1. j> i/dy 
«+( vu+i, j i+vn + u j-n-vi i » j)-vi i » j-i ) ) /dx) 

GO TO 451 

0453 IK FlI + l, J I .ME .4) GO TO 454 

IF(FII.JfI). FTi « 4 I GO TO 451 

PHI 1 1, J )=-.5*FTA( T , J )«( lUI I ,.1»1)+U( 1-1 . J+l) -Ul t . J 1-U( 1-1 . JU /OY 

*<■( VI 1, JHVI1 «J-1 I -VII — 1 . Jl-Vll-1, J— 1 ) 1/PX ) 

GO TO 451 

04 54 I F< F ( 1 , J+l ) . E 1 .4 > GO TO 451 

PHI ( I . 11 = 2.0*5141 I . J l+l V( 1 ,J)-V( I ,J- 1) )/DY 

• Gil TO 451 

045? 1HF1IM,J).W=. 4) Gl) TC 1 4 55 

if (Mi .j+ii .r.'' : .4i on to 455 

1F|F( 1-1. 11 .1 ).4I GO TO 451 

PUT ( 1 , JHO. Air,J)*UUlI»J)«-N(t-l t J)-U(l,J-l)-U(l-l.J-lll/OY 

» + l VI 1, MfVl I, J-j ) -V ( 1 - 1 , J I -VI l-l, J-H ) /OX) 

GO TO 4 51 

0-.56 IMF 1 I -1 » J) .1:0 .4 ) GO Tt) 4 51 

PH Ml, J) = 2.0<-LrAI 1 , jJJgjUU .J )-U( 1-1 , J I l/OX 



£0 . Jt LA.? I 

045b I F I F 1 I.J + l 1.NF.4) GO TP 45 7 

Ji,lT' U- 1 1 J 1 .j. }.4 ) J'>0_Jii_A5\ 

PHlT l. j ) = P.ntLTAH .~J ) 0 (V! f.,1 I -VI 1 . J-l) )/))Y 

0 ) T ■ 1 ?■ j 1 

04 $p ph i ( l , j! r a < Li J im ( :ii i,j >+u i L:L.JH)I i,j-i )-utl- 1 . J-i)i/py 

»«-! V(J_tl ,_J ) + V ( I *■ 1 , J - 1 ) -V ( I.J)-V(1.J-1 l l/r xt 

£fi_T ' l ± 5 l 

0437 IFIFI [-l. Jl .fii-. A) GU T< I 451 

. P'jH T, M = ?.0*l : TAM,J]«lU(I.JHJ(I-l,;))/DX 

0431 CON T [ H JF 


£ ST04F VAPTArtIFS FUR RF— FVC I ING 


naiTF116l'J-.\/. c rA,PHt.P 

on tr 403 3 


C Rf. -ESTABLISH VARIABLES FCR Rif— CYCL P 


A 036 REA0116)U.V.ETA,PHI,R 

1FIFQF, 16)4036,4036 

4018 P£v,'[vQ l r> 


C CALCULATE EST [ FAT E OF FINAL WHEEL VELOCITY 


UG=Uv-i-H0XFl SFXG/.P'l 1*0T 
VC=VH-HGY + ( STYG/ I* ) ) *0T 
fx={ i (in;tUK)/7. hot ) /nx 

FY-( { ( jr.+vu)/?. )»11T) /nv 


C CELL VARIABLES AFTER SUFTt NG COPPIHATES 


!lai 

N^l 

TF{ FX.GE.0.0 ) GO TO 417 

M»-l 

0417 tFIFY.Gr.O.O) GO TO AIR 

N = — 1 

0410 FX= AGSt FX 1 

F Y= ,ViS< FY > 

m ‘t?') J- ? i JftARl 

DO 4 70 1=2-1 GAk I 

1 F1F1 1, 11 ,NF. j.A-in.Ft 1. ) 1 .NF. j 1 GU TO 420 

[FIFE (I. JUNE. 2 1 00 TP 421 

IFIFU-I-I, J). !?•;.)] GU TO 422 

PHHltl,! l-PFI 1 i , JH-GX»[)X<-7. *ETA( 1 1 J 1 *U ( 1 - 1 , J 1 /O X 
90 TO 423 

0422 IFt F[ I + L. J) .M-.5 1 GO TO 423 

PHI 1 I tl, J)=PHIU , Jl-lSFYG/wr )*DX + 2.*ETA( 1 , J 1 <- 1 U ( 1-1 , J l-UWl/OX 

0423 1F(F( 1-1 , J) .NE.l) GO Til 424 

PHI 1 1—1, J I=?M1I It J}— GX*DX-2.* C TA( I, J) *U< I, J1 /OX 

GO TP 425 

0424 I F( F( 1-1 , J 1 .NE.51 GQ TO 429 

PHI ( 1-1 , J ) = PHI (I, jl4<SFXG/WM)*DX42.AFTAl I, J HMUH-Ul I ■ J ) ) /OX 

90 TO 429 ; 

0420 IF! F( 1-1 , Jl .FiE.6 1 GO TO 425 

PHI! I-1,J)=PHI1 I. J)~GX*OX 

0425 IFIFCI, J + l) .ME. 11 GO TO 426 

PH1M.J+1 [=PH! ( i , J H-GYFUY+2.*ETA1 I . J)*V( I, J-l ) /DY 

GO TP 42 7 

0426 IFIFI I, J»l) .ME. 5] Go TO 427 

PHI t I , J + l 1=°H1 I I. J1-1SFYG/HF- )*t)V + ?.*STAI I, J 1*1 VI I , J- 1 )-VH 1 /O Y 

0427 If! FI I.J-l).Nc.H 90 TO 428 

PH Id , J-l )=PHI ( t , J)~GY*UY-2.=>ETAI I. J1*V( I,J 1 /DY 

GO TO 421 

04 2 S 1 F 1 F ( 1 1 I- 1 1 -lie . ? 1 GO TO 421 

PHI ( I . J-l I -PH I I I ■ J 1 +1 SFYC/rt* )»t>Y + 2.«ETA( I,Jl*(VF-V(IiJ>) /DY 

0421 PHI 1 I ,-)) = ! t .-FX )«l l.-FY )*PHl < [, J) *FX*( 1.-FY/2. HPHI 1 MM, J) 
s+FYd 1.-FX/2, 1*PH1 I I, J+H ) 

0420 CONTI NiJE 


C REGION 40 °HES SURE DENSITY CALCULAT TOMS 


400 tr.=o 

40 1 CON^AX-O.O 

up 402 I TFH = 1 .10 

no 403 J=2,JRAK1 

liil 403 I=2,IP.A4I 

!HH I, II .NF.21 10 T 2 <*03 

I F I r- r 1 1 . )~l . 1-iV . > 1 GU TP 405 

IFIFIH-l.JI.Ne.il PQ in 407 

PHIU + HJHPHKI. JH-GX»I'X»2.»FTA( I,J1*IJI i-l,j)/nx 

90 m 4 06 

0407 I F(1 < IU.J 1 .NF. 51 90 TP -.06 

P ut ( I +1, 1 > — PH III. J)-tSFXG/Wl«PX + 2.*FTAt 1 1 J) MUI 1-1, Jl-UW) /PX 

04 16 trir n-i..n .m=. i ) go in 409 



PHI ! Lzl» J ) nG>S*£ifcZ*J? LUUJ i £ J .IUJJL j >UZPX_ 

I GO. fn 40>i 

040G i'hi'iT-i.Tii.m ? _5 ) 4.15 

ElU! Lz t . J.) l£U.LLbJJ_+_LSLx >VW ) «n x I , J )«l UW-'JH , J) t/»X 

on 11 40 •< 

0415 IF U ( I_-l • J on Til 403 

pH i('i-7t j 1 =7m r7~t . j i-gx^o x 

04 )ri 1 p 1 r 1 i* )nj o_jri 4jj 

._e.HL{_iZu- u=p'r?r tJj.jr»i^j '>y* z,*£iA.ii.. j..i "V-U-..J-1 > /fli 

on in 410 

04 1 1 1 PI I I I. H-Il .KC .4) on TP 410 1 

f hi ( 1 . j« 1 ) = ph 1 1 1 . j ) -is f rr./wn nnY+J .* £ Tfli 1 . J )»t vi 1 . j-l i-vwi/oy 

0 4 1, 0 -JLE.LE ( JuTt-LL , m , J-LJXU_U1_VL3 I 

PH I t I ■ ■)— l ) = PH I I I. J]— GY»l'Y— <;■ »f-TAI 1 »j i*vt t . J i/iiy 

on TV 40b ~ 

0413 IFird . J-l). ^.51 on TO 40 0 

phi 1 i , 1—1 i = phi ( 1 , .1 i+isfyg/w 1 )<-hy<-7.*{:t<u n j )<m v^-v i 1 . ji ) /nv 

4 05 PHI is; /I = { PHI I1H . JMPH1 I 1-1 ,JI l/IMIXS) 

<■+ 1 phi n . j 1 1 h-phi ( i . j-n ) /i /ftOYsi+p u . j i n 

' ini 1FN. G E.101 0.1 TO 412 

COK'I J=4flS (PH inert- PH I t I . J 1 I n APS! PHI HF'J 1 + 4 OS I PHI C 1. J) ) 

*4-0. ?5< (III I , J 1 Hit 1 - I , J ) ) 3 40. Pgff (V 1 I , jTfV I I , J — 7 ) 1 S&24GVYT4QXXR1 

' IFICXIJ-LT.CUNHAXI OH Tf 4 12 

i CTiim^jnij 

■ 41? PHI < r. 11 = PH1i- , FJ 

403 awn mp 

4 02 lf.= If»l 

if t \.no-i r . pi on rn 414 

IF! IC.P3.Mir)->0Q.40t 

0414 iricivr-nx.GT.r.iirivi on to 401 

c 

C PcGin-I 50 VELOCITY CALCULATION 

C’ 

05-JO M=fl 

• I CPlJNTaO 

: 00 511 J-= 2. JR AR 1 

’ no 511 l-g.in&pl 

N=N+1 

IFIFI I, J) .NF.51 00 TC 503 

TU(X )=!I0 

T V t M = VG 

PHI(l.Jl=0. - 

GO TO 53? 

0503 CONTI MUg 

IFIFI I, .11 .NE.;.AXn.FII,J).NL.3) GO TO 520 

UAflOVE-U( 1 1 J-l 1 

U-JELnu-UI 1 ,J+l) 

vt^vi i , j-n 

ve=v< i,j) 

' VTR = V( 1 + 1 , J-I 1 

VSft = Vt I»l. J 1 

IFIFI rn, J) » EO .4 ) GO TO 512 

IFIFI H-1.J1.E-T.il GO TO 513 

IFIFI H-l.JI.FO.5) on TO 501 

I H FI 1+1, J4-1 ) , c 0. l.OP. FI I, J4-1 ). e 0.1) GO TO 514 

IFIFI 1 + 1. J+l) .EQ.5.PP.F 1 I . J + l I.E0.51 GOTO 5140 

IFIFI H-1,J + 1).NF.4.1jK.FI I ■ J+ H .NS . 4 I CP TO 516 

U3ELiH=')( I .J> 

_J GO TO 5lo 

- 512 TUI XI =■!! I . J)+GV<-HT 

GO TO 51 T 

0513 TUIMleJ. 

GO TP 517 

0501 TUIN1=JG 

GO TO 517 

. 514 U'*FLPh=-U< I ,J) 

I V3 = 0. 

VHfi = 0. . 

GO TO 516 ; 

5140 UJFLCHs 3. SUW-UI ItJ) 

vn=vw 

ypfr=vrt 

0516 IF! FI 1 + 1 , J- 1 > ,~0.1 .OR.FI I. J-l I .EQ. 1) GOTO 518 

IFIFI 1»1 «J-tl.gQ.5.nx.F( 1, J-l) ,t:0. 5> GOTO 5180 

1 lnrim . j-i ) .NL-.4«t)fi.F< 1 . j-ii.ne.4I go to s?o 

UAflOVE=J( 1 . J) 

GO I 1-1 520 

518 UAH )VF=-U( I , J 1 

vr=o. 

vrr=o. 

GO TO 5?0 

5100 OAbflVf = ?.«U-V-U( I, J1 

VT»V W 

VT.l - V W 

SdO SUM I-.IAII.JKIIJIHI.JIH; I I - 1 , J 1 -3 .n»m I, J 1 ]/f )X S» 

+ IUP E 1 nil+UA GOVE -? . 0 F II 1 1 . J 1 1 /DY S 1 +1 PHI I 1 .J1-PH1 ( 1 + 1 , J1 1 /11X + GX 


Oi&A 


U% 

Vt 

'tt 


CD 


Q. 


vc) ^ 
a ^ 
o. 
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SLIM 2-1 iMrtflVr + UI J« Jl > * < VTH + VII /4. -lUKFt f’W VI ) ( 1 , Jl ) * < VEP + V» )/4 . 

'57 Y~ Yu ( nT = ji i , j) +*ji*( siimj iijii, j j +u< i- 1 . J >)*» 2 /nx' 

»-n ./F»un 1 + 1, !)+’)( I , J ) ) * +2/u x ) 

AUL^uVC^.I Lri+JJ I 

{■'111 *W I t A. I 11 


— 


LM “ 

■ 1 ~ v ^ i r t i.'i.x 

u( l ,J> 







U:W' 

-111 I.J+L) 







ui = 

U(I-l.J) 







in I, 

*uu-i . j+n .. 







_LF ( 

I: 1 [,J+1).F').4I r.0 TO 571 





iKfl i, .i+i i.m.n r.r re 5?? 

IFIH 1 i 1 + 1 1 .F’J.5) CP TO 507 



IFIFI 1 +1 ..1+1 l.rO. 1 .UP.) t 1 + I..II.FM. 1 ) 

_C£L 

rn 

523 




1FIFI I + l.J + l 1.F0.5.0P.M 1 + 1, Jl .F0.5) 

GO 

TO 

5230 




IFIFI I +1, J + I ).NF.4.0F.F( 1 + 1. J) .NF.4 ) 

nn 

TO 

524 


vi<ight=v< i , j i 



*0 

T° >?\ 






5210 

V>IGHT=7.#VW-VlI. Jl 







m- 

UW 







U8F 

=IIW 







g n 

TO 524 



* 



■ 521 

rviNi-vi i . ji +gy*ot 







CO 

TO 53? 






052? 

TV 1 PC ) - 0. 







GU 

TO 532 






0502 

TVIN) = VG 







GO 

Tr) 532 






_ 5?2_ 

VH ICHT=-VC l« J) 







1)4 = 

3. 







IJ!}5 

= 0. 






0524 

IF (Fll-I.j+U.eo.l. OR. FI [-l.Jl.FO.il 

GO 

TO 

525 




IFIFI I-1..I + ! I.F0.5.UR.FM-1.JI .F0.5) 

on 

TD 

5750 


IF( F( 1-1 . J ) . F 0,6 1 CO TO 507 



IFIFI I-l.J+1 l.\JF.4.PP.FI 1-1.JI.NE.4) 

GO 

TO 

526 


VLFFT=V( I , J ] 



Gil 

TC 526 






52 50 

vlfft=>.fvw-vi r ..i) 







UL = 

UW 







U*L 

=uw 







GO 

TJ 576 






525 

_V.L\<rFT»-V( I . Jl 







IIL = 

0* 







U3L 

=0. 





- 


GO 

TO 326 




- 


0507 

VLCFT^Vf I , J) 







UL = 

0. 







U3L 

=0. 






526 

SU M 1 =FT 4 1 I ,.!]*< ( VI 1 , J + l )+V 1 1 , J-ll-2. 

o*vi i ,j n / dy s 



* + l VK JUrlT+VLFK-2 . J*V( 1, J 1 I/DXSI+IFHI ( 1 ,J1 - 

“HI 1 I i J+l) l/DY+GY 


SUX2 = (:JL+'J5L 1* (VI I , Jl +VL cFT 1 /4 (Ufi+USC) * ( V 1 I . J) +VP IGHT1/4. 


528 

TV(N> = VI I , J ) +OT<=l SUM l+SU’^P/OX+O. 25* ( VI I 

, J l+VI I ,J-1))#*Z/DY 




F-0. 

2 5*1 VI I ..II +VI I . J + l 1 1**7/0Y1 







GO 

TO S3? 






0529 

IFIFII+1.JJ.N5.5) GO TO 504 





TUINI=UG 



GO 

rn 






0504 

TUI N ) =U ( I , J ) 






0505 

Irl 

F1I.J+11.ME.5) GO TO 506 





TV(\')=VG 



GO 

TO 532 






05 06 

TV(M = VI[,J) 






53? 

IFIN.l T . 1 020 1 GO TO 511 





UP I TF ( 1 l ! ( TUI H> , TV(N) ,15= 1, 1020) 

ICeUNTKCOUNT + J 

N=n 


51 1 

CON T 1( V UF 





I T* C H } IT HN) , TV(M ) ,N=1 ,NIVI J ) 

7EWIMD 11 



CO 

TO 540 






540 

N = 0 






*4=0 



DO 

541 1 = 2, JM"! 







no 

541 I=2.1B4Fl 





N-N + 1 

iriN.HO. I ) CO TO 542 


544 

UI I 

. J 1 = T U 1 N 1 







VI [ 

, j) *rvim 







GM 

TP ‘5 <}<1 






5 

IK 

(.09. I COUNT 1 GU TO 54 3 







» r A 

0(11) 1TIJINI , TV(N) ,N=t. 10701 





M-MM 



GO 

T’T S44 






54.1 

°eAoi in ( turn , TviNi ,n=i ,mm r j> 







rn 

rn Wk 







IK I -l.l T. 1020) GO TO 541 





N*0 
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— 5.7.1 XO£)lIN<i£ 

nil 574 J = 1 i J8 A » 

U U .• J ].£9.f I 

1 EiJtS.tiJt.N5j.U_ TO 5 75 

V(1.J)-V(2,J1 

r, o 1 1 ; fo 

-P.5 ,J\.::rV..U,JA 

-aJL7i>„UU_‘) V.'jJ ) zo., 

5m_.L.f!Nm LlL 

on t ■ kkft 

060*; CONTINUE 


r COMPUTE FORCES BETWEEN GHEFL AND SOIL 


sfx=o. ; 

$rv-n. 

ASSIGN 653 Id KRET 

on 6 S 1 j=2,jeAui 

on 651 l^lj J JAKJ 

1FIFH t, J) .NE.2) GO in o51 

lF(F(I+l.J].Nr.S.ANn.H I-i , J] .ME. 5. ANO.Fl I , J+l I . NE . 5 ■ 

*A?in.F( I. J-l) -NF.5I go TO 651 

SAVS=ETA1 1 . J ) 

KFr-p-rc I [ - J ) 

IU F( t.J) .NF.? .AMI). Ft I. J) .Nt .3) 651,611 

06 5 3 SU , <1 = ?.*0Y»FTAI I . i ) <- H) i I , J ) -U (I - 1 . J ) 1 /nx 

1FI S'J'U.LT.O.-)) 00 TO 654 

5 UHf_n, 

0 6 54 Slin?= ■’.FQX»FTA( i . J ) *( V { 1 . J I -V (J . J- 1 ) )/0Y 

IF I S'H?. 11.0.01 i;n TP 655 

SUM 2=0 » 

0655 I FI PHI ( 1 ■ J) .LT.0,0 ) CP TO 656 

Siil l 3 = PHI( I,J) 4=010 S*DX 

SU.-t4-°Ml ( I , J)^KHOS*0Y 

GO TO a57 

0656 SUN1=0. 

SUP4=0. 

0657 SU;-'5=TAUXY( 1 , J)*')X 

SU.’i6=t AOXY { [. JU-QY 

FT AH . J) =sayg 

FC ( 1 . J) = K££P 

■ IF( F( I . Jf I ) . NE.5 ) GO TO 660 

SFX-SFX-SUM5 

SFY=SFY-S'J'12 + 3>m 

0560 I F ( F I I » J- 1 ) . N- . 5 ) GO TO 658 

SFX=SFX+SUK5 

SFY-SFY+SU’^-SOMl 

0658 IF(FtU-l,J).A£.5) GO TO 659 

SFX^SFX-SUI'1 + SUY4 

STY-SFY-SUYS 

0659 IF(F( 1-1 . J) . ‘Ic . 5 ) GO TO 651 

SFX=SEX»SW1-Sin» 

.SFY = SFY»SIH6 

0651 CCNTl'i-lc 

if i ifsym.he. n go to 66i 

SFX-O. 

SFY^P.FSF Y 

06a I CONTINUE 


C RFGU1N 80A- CHECK t ST IF* AT F OF rtF'EEL VFlOCITY AGAINST 

C CALCULATED VAHJF, ADJUST ESTIMATE AriO RF-CYCLE IF NFCESSARY 

C CCNPQIF MEN hhFEC VFLnCITV BY I FiPULSE — NOHFNTUK PKICtPLE 


UWN-UWHGXtSF</W’M *QT 

VwN=VK*-( 0 Y*SF Y /a'N I X-QT C 

1F( IFSYM.EO.l ] GO T£ 803 

IFI SFX.\'E .0.0) GO TO 810 

[F(SFXG. c O.O.iU QO TO 822 _ 

NORgQ 

SFXC'Q. ~ 

SFYG=0. 

GO TP VI36 

032? ERFflRU = 0. 

'GO TO 501 

0810 FRRO'iH = ( [SFX-SFXr,)/SFXKIOO. 

IF| A6S tiRRORU) . I T .0.0001 I GO TO 600 

IHNjK.GE.ll GO TO 805 

0800 SrxiH.PX 

SFXOl-SFXO 

SFXOH SFX-fSFXGI/2 . 

G'l TO 10 i 

0805 XI = SFXG ~ 

SI'XG- ( I SFXl*SrXG)-l SFXYSFXGI) I YISFXG-SFXr,I»SFXl~SFX ) 

SFX l-SFX 

SFXG~l=Xt 
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08 03 I F I SKY .NF . 0.0 I C.U TO _8 1.0 

I f ISFYG.EO.Q.i)} i}5 TO 8ZS 

i\ll]K = ff ~ 

sfxg »o. . _ . 

S£Y Q--Q , 

Ca ) T i ) hOJp ____ 

08 38 

O U TO .807 

081') l.' < UllPV-l i~s F Y-SI- Y'~, 1/81 Y)- »100. 

[riA'vs (epgrpy L. i r.o.o noi ) fill to p .06 

li rnr.OF.i ) f.n rn .ioh 

nnp( .SFY1=8IY 

SFVGi = SfYr, 

SfYG=l SFY»SfVG)/P. | 

c.n rg no 7 

. .iIita2-Vl^.S.D23 : 

SI-YO=HSrYlMSrYO)-( SFYSSFYC,!) l/(SFYG-SFYr.l*SFYl-SFY) 

SFYt=SFY 

SFYG1-Y1 

0307 IF ( AO SCM U' RU) .LT. 0.0001. AND. ABSIERR0RV1 .IT. 0.0001) CO TO 804 

IFI ABSt IXP-T ) .GT.I OT/Z. ) > GO Tt) 809 

V(H I IF (6,14 I T -I.MN .VMM .UP, VC. S F X , SF Y . c RRORI I ■ ERRUkV . IC 

W° IT< ( »>. ? I ISFXG. SF YG 

0808 T F { l\jQP .08.10) GC~ Til 804 

MO 8 =0-01 +1 •_ 

IfmfH.NF.3) oo TO 4036 

N 1 C=IC 

GCI T~i 8036 

0 304 CI1NT1 MUF 


£ REGION 30- CALCULATE MOVEMENT OF WHEEL AND TEST FGR SHIFT 

T. GREATER THAU ALLOWED 

X RE-CYCLE WITH REDUCED OT IF NECESSARY 


OXSF=( (lHtl)G)/7.l*0T 

0VSF = ( ( VW+VG) n. ) ^r.T 

I S F= ( 4 .X'OXSF ) / OX 

JSF = (4.»ilYSF)/t>Y 

IFdSF.EQ.O.AVO.JSF.EO.O} GO TO 857 

W8ITEIS.17IT 

TgT-OT/4. 

tn = .7 5*r>T 

DTCP=.75<=DTCP 

0TPP=.75*0TPP 

0TSP=.75*flTSP 

TCP=T 

TPP^r 

TSP-T 

.00 R = 0 

SFXG-O. 

SFY0=0. 

M=UTF(6tL6)DTtPTCPtOTPP.DTSP 

GO TO 4036 

857 CCSMTI»’U£ 

IF ( UW.EO.O.Q] GC TO 6870 

ALPHAIJ=UG/im 

6870 1FIVW.SQ.0.0) GC TO 6888 

ALPHAV = VC/VH 

6368 CDNTIF Ur 

UH=UG 

VW= VG 

SFX0=SFX 

SFYG=S FY 

SOX3F = SOXSF*OXSF 

SOYSF-SOYSF-HIYSF 


C REGION 60 COMPUTE STRESS TFNSQR OF FACH CELL 

C Fmo PRIATtPAL STRFSSES jM> DIRECTION 

C TEST YIELD CRITERIA 


MU- 1 

STNMAX^O. 

ASSIGN 615 TO XRST 

0600 OH 601 J»2, J8AR1 

on LOT t*2,I8A3L 

GO TO 16 11,612,6131 KM C 

0611 FCI I. J I =0 

tau<yi i , j ) - o . o r 

IF! FI UJ I .Nr. .7 .A MD.F1 I,J ) .NE.3) GO TO 601 

U5P-C.K I , J + 1 n-tn I , J) 

npL - 1 n ( 1-1 . JH ) »IJ( l-l.JII/Z, 

UA 9 *CU 1 , J ) Hll I . J-l I )/Z. ~ 

ual- iu i i-i, j ) m i i-i,j-m/z. 

VI Ad VI 1 ,J-|)>VU-I,J-H l/Z. 

VLO= (V 1 1 , -J 1 »Vl 1-l. JH/ P . 

VRA-I VI I «■! , J-l I + VII , J-F) M2. 
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. 0:1 . rn_ 5 Q 3 

0£ Qjr t £.{ t- Lb J *■ I ) . n f . a l _r, o TCI 60 9 

1 Lt 1 FU - b J i .1 ~i , -t L 

6UlU (0 lb J) iJ-ll=mJL--b,bJJJ/n 

!>UM>=lVPft-V(. A>V'M-VLHI/t ? .»I)X ) 

OH III 6 0? 

QMO lz.l.J.-in /0y 

I V LLM I J!_*_Y U±.bJ=l !-.V.tJb-JJb;YXJbJ-lJJ m 

r.n rn 6 n i 

06Q<? i fif i i-i . .n.ivt - . 'ii gq rn ftm 

S'jr'1-i'.nw-nAo HJUL-II*! )/( ? .»QV1 

Siiil 6=.>*< v ( bb-Jbvi i»bJ-l >-vu . J)_.v( b J-l U./OX 

GiL £ii_sLU 

Q6jr sumi = i n . r-'.iap»‘ 1 pi — u a i i /t ?.«nvi 

sux v- 1 yp. a-vYa+ v h»-v "l h )/12 .»nxi 

06.03 SU.jJ=-2..J?(.E XAJ i.i J 1**Z 1£X U <U(I x J I zll i L=. l .JH/P X)**2» 

*H (IVI I . Jl-V( It J-l n/DY) **^) + .25*< ( SU"H-SUV)^n 

STPN£O.JL?( SU “1 + S'J M2) 

STRAIN ( I I = /.»S ( S TRN I 

I F ( ST *4 IK ' t 1 ) .1T.STNMAX) fin TO 629 

S TN‘«MX-STRAI Ml ) 

nispq rmTiPOF 

TTA(.i , .n° aAg 

TAUXYt l,JU?.9*FTA<I,.))gSTKN 

EH T . J > = 1 

r.n to KR r t - 

615 FK f 1 y.M=FK1 

IX.' t s un. LI . (FjUI. J ) *■> .? U01 J T Q 10.1 _ 

IFtSTPAl M 1 ) . Lfcj.OAyAti ) GILT 1 ! OQI 

rail <PI INC flMPTS. MMCVS.M^AX. AS. XiCVLMS.STKAN.STRFS. STRAIN. FPS, 

*PROX! Y. STRESS. SSI .SS2.S2.S3. I1ELY.H. IP. IT01 

rTftll, J1^0.i*‘STRrSSn, n/STRAIUH) 

IFtSTR'I.GF.O.O I r.n Tfi nrtO 

Tflnxm .Jb -S Jft. p$ 5 tl .«.U 

r. 1 TO -.31 

6R0 TAUXV1 I .JI=STRFSS( 1. 11 

6.01- t P=,LbJ 

FC ( 1 . J 1 » 1 

on to 50i 

Obi? TFI FI I . J) .\F.7. AND.M T. I ) .ME. 3) RH TO 675 

SIG-'AXt t,J)=?. Jt>TA(I.J)*(U( UI-1HI-I.JI )/DX 

SO TH 601 

0675 ETAH..tl = 0. 


STGRAXt I , J1=0.0 

GO TO bOl 

0613 If-IFI r,J).Ne.2.AMi).F{I,J) .MF.3) Gil TO 677 

S1GKAYI I . J)°?.7*PTAt I..H*(V{ I, J)-V< I, J-ll 1/PV 

r,n to 601 

0677 FT A ( 1. J) = 0. 

S 1 G M AY ( 1 , J1=0.0 

601 CONTINUE 


£ STORE STRESSES ON TAPF 


IFtA35tTSP-Tl.GT.tCT/2-.) 1 GO TO 673 

OP TfM670.67T.677) wy 

06 70 UR I TFI 9 IT 

WRITF191TAUXY 

GO TO 673 

0671 WR 1 TEtOlSlGOAX 

GO TQ 67) 

0672 WRITF19I5IG 1AY 

FMO FILE S. 

TSP=TSP+DTSP 

NSP=WSP«- 1 

WP TTg(6,251T.NSP 

0673 C'JNTl -lU^ 

mm=m‘h-i 

IF(r!M.LE.3) GO TO 600 


C CAICULATE THE CFII C ISCREPAMC1 ES 


I 0 U 6 X = 3 

JUPAX-O 

DMA X = 0. 

fl'.l 6036 J = 2,JRAP1 

60 6015 1-2.IJA F1 

(U I . J >=0 

At)PAX*0. 

tFIFI I.JI.6F.P.A40.FU.J1.MF.»1 CO TO 6035 

ft i .ji=((iii i . n-ut i-i , ji i/nxi + ttvt i ,j)-v< r, j-i n/QY) 

AOMAMOt i , J 1 

1 F| 1 .( 1 ,.l) .1 r. o .01 AIV 1 AX=-Lt t ,J) 

T F t o«<x.r.F. AiriAX) go tii 4035 

IIMAX = AtHAX 

t OM AX - 1 



jro* a Xfj 

401s cow n Ml IE 

X. - „ „ ,_ 

C. r.L!I_£i!tL.XLM L-XL)-.tlEi JO_i: f. LLVXBi APJXS_U_N JUM. 

X 

if (.aos { icp-jko r . iqt/?.u go rn 7?o 

KiiJJ.L.lijj.23 ) 

U£LAO^_Jfi'Ui jJUF. 

.BIT AlUi— UKJLS-iHli 


ysi.q . -HA .x.t ,|, iT ,j[ 1 , J J IH 'x yi 1 , Jl i. H 1 r JJ j JE Ui J ).i. K U -. 

40 33 CtlMTKOi: 

rrin.^MlT 

07 20 COUTIMUE 


Pt-OTON 70 PAUTtri F MI1V6MHXT 


NN» 1070 


fin 721 H=1,HPARL 


IF t M, LE.OPAR) GO TO 722 
WjgMPRR 


07 7? REAP 1 14) I TXKINI . TYK( N 1 .N = l .NN> 


no /n I -M = 1 , 1 ’\' 


xmm = T«iri) 


YK I M =7YKIM 


I=XK(M)/nx*2.0 


J=YK1M)/0Y«-2.Q 

JE 1= l 

JLLrJ 


FX^XXIN) /DX*?.0-FT 

FY=YKIM)/DY*2.0-FJ 

|fini,J).MI..?.Mn.F(l..H.Mg.1) no TO 701 

S<glX(fl-Xit(M t/i)X 

-LU^MItJ.. J.)... _____ — 

,» ? = UUt J > 

iFiry.i t.o. 9) on to 703 

SV=|VPL( J)-YKINI ) /DY 
i-i, )»ti 
U4 ~U I 1,JH) 

I F I F c ( T-l . J) .PE- I I GO TO 704 

m=un 

U 1 - u w 

r , 0 TO 705 


0704 IFIFEU-ltJ+H.NE.I) SO TO 705 

LI3-U1 

0705 IF(FEI H-l.J).>lE.l) GO TO 706 

U 4=017 

U 2 ~ 1 J W 

GO TO 707 

070g IF(FE(M-UJtI).NF.I) GO TO 707 

t)4 = U2 

Gn TO 707 

0703 SY=(YK(N1-Y3U J-I ) )/0Y 

IJ3=U( I — 1-J-l) 

U4 = U(1,J-1) 

IFIF2I l-I-.H .NE. ] ) CO TO 70S 

U3 = LM 

u l = uv> 

GO TO 704 

070S'~IF!Fr[ l-l.j-l l.tl-F.i 1 GO Tfl 709 

U3=Ul 

0709 IFI FEl Itl.J I.NF. I I GO TO 710 

IJ4 = Ur) 

02= UK 
on T1 707 

0710 1FTFFI I + l.J-ll.NE.l) GO TO 707 

1)4 = L)2 

707 tK=(0.5*SX)»(0.5 + SY)<=Ut» <n.5-SXItt|0.5fSV)*U2 
*-t-(0.5 + SX)*IO. j-SY )*U3+|6.5-SX)*(0.5-SY )»04 

X K(l'.'l = X i<IN)»UK<'DT 

SY=(Yt Jl-YMMU/OY 

Vl-VI I . J- t ) 

V2 = V I 1 ,J1 

TF1FX.LT. 0.9) on rr 711 

SX = ( X2 LI I)-XK1-V))/DX 

V 3 =V ( I + 1, J-n 

V4=VU + I. J) 

IFCFFj I, J-n .NEi. t I GO TO 712 

V <= VW 

VI = VW 

Gil r< 713 

07f2~ IITfiTk-iTj-I I.N g .l) GO TO 713 

V 3= V 1 

0713 IFIFFI T.jni.MC.l > 0,0 TO 714 

V4=_vu 

V2--VW 
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on r n ns 

0 - 71.4 tf[ f H. jiM.,j»l.l.Nf ,i)_t;o.rn 71 s 

V 4 _ 5 .v: ! 

r,n to 7 15 

...Zll . 1 1 ZD.X 

V 3 - V I 1 - 1 J.-IJ 

Vir tf'< I - 1 . J 1 

If (FP( l.J-l l.NE.l) GO TO 716 

Vl = tfK 

V l = V il 

i;n T'l 717 

mu, iFirt ( i-t, j-ii.ne. n go to 717 

su=.YJ 

0717 fFtFFlt.JHl.NE.l) BO TO 710 

V 4 = VW 

S/?*JH 

on to 715 

0710 IFtFHI-l.J+lKNE.ll GO TO 715 

V 4 = V? 

715 l'K= 10 . 5 »S X) Sl 0 . 54 -SY)aVl<-I 0.5 + SXl*( 0 .. 5 - 5 y)<‘V 2 
ft+f 0. 5 -SXI » ( 0 . 54 - 5 Y) 3 V 3 +( 0 .S-SX)« 10 . 5 -SY)«V 4 
vk i m = vk r.MittfKanr 

701 COUTH- .jf 

HOT TM 15 I IXKim,YK(H).N=l.NM) 

7 ? 1 GONTINUF 

on rrn 14 

3 EWI.T 0 15 


X COMPUTE COnPOT HATES OF PARTICLES AFTER SHIFTING 


on 047 i»i,jP4s ; 

00 34.7 I = 1,1 BAR 

KC C 1. J )-0 

03-4 7 CONTINUE 

KPUM =0 

KPL P-0 i 

_ tufcoaza 

00 e4l H=1 ,'IPARl 

1F1M.LE.RPAP I GO TO 342 

1JN = PPPR 

0347 3 EA 0 I I 3 1 1 XK 1 *3 1 , YK 1 N I « M— 1 . NH > 

no an n=i.w 

1F( (Y<<(.4)-DXSPl.LT.O.O.nR.(VKfW)-OYSF).CT-YT) GC -TO 811 

KPLU«'=KPLUH-H - 

KPLP=KPLM-fl 

TXK (KPLUP )=XK(N)-PXSF i_ 

TYK( KPLUtt)=YK(H)-l>YSF 

1 = TX K 1 KPLUP ) /0X+?,0 

J=T T 1 K t KPL UP ) /PY+2.0 

IF1F11.J) ♦ll c .5l 00 TO 344 ’ 

n^dXKlKPLU.PlfPXSFl/OXf?, 

IF1F1 n.J l.CfJ.51 GO TO 345 

TXK 1 KPLUM) =TXK( KP LU. V ) 4PX SF 

1 = 1 1 

GO TO 344 

084 S J1=1TYKIKPL I JP|4-0VSF1 7PY + 2, 

1F1 F(1,JJ) .1:0.51 GO TC 046 

~ TYK~lKpl~UMl~=TYKiKPLX)>'l-fDVSF~' ~ ~ ~ ~ ~ ' 

J=J.) 

GO TO 344 

0846 HP 1 TE 1 6,60) 

H3I TF(o t 5 9)1. Ji TXKlKPLWdTVKIKPLUMltFU, J) ,FE( t,Jl 

GO TC 0Q1 

0644 CONTINUE 

KC< I,J)=~KC( I.JIH 

tFlKPLH.E0.U201 GO TO 312 

go t.t m ; 

812 HRIHUA) (TXK1N) .TYK1N) ,N=1 ,10201 

KPLUM = Q 

811 CONTK'UF 

841 C0NT1M.1E 

813 N = K PLUM ■ 


K-K PLM 

C CCM8TL THE OlSPLaf-FMENT OF SHIFTING COPKOlNArES LESS THAN 
C HALF CELL HI MTU EACH TIKF 


CiNX = \X 

CHV = >IV 

XKMx-’n^icT’x-i.oH'nxK 

V KMY=VJHCWY-1.0)--PYK 

S IX-XP-f-GOXSF 

SMY = YT«-SOYS F 

M=< S<X-XKMX1 /n.XK-TSFPAR 

1F1H) 326,0 10 ,676 

0328 XSF PAP.XSFPAP— 1 . 
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Q0 _T:l_d L6 : 

0876 f.UNT I MUE ; 

v p = vo -in y §f ; 

xp*xni-ir .iK+xsrp ak ) adxk-soxsf 

no rii rry=i.MY ... 

N^M+_1 __ 

T XKtM) = <P ; 

JVK ( ii ) »YP 

illxKCi!Y/r)*tZ,.0 

ihjlyk U'UV.DYt^j 

KC( t. j )-<r.( i « j m 

ifim.l r . 10201 go ro ai5 

wa ne U 4) (txk(m),tyk(W)tN=i,io?o) 

N -0 

815 K-K ■*- 1 

yp = Yi»*riYK ; 

814 COMTTwjjF 

XSFPAR=XSFPAP*1 .1 

016 NN=1 SmY-VKNY I /OYK-YSFPAR ____ 

IF ( N’l) 321.8000.877 

‘ 0839 YSFPAF -YSFPAR— L. ’ 

Gil T-,1 8000 1 

00/7 CGAT1MJF : 

YP = Y0+1CMY<-VSFPA8) A PYK~S DY SF 

xp-xo-KSFPfin^nxK-snxsp 

HlV 818 ICX=1 .MX 

M=N+1 

TXK(H) = XP 

TYK 1M)=VP 

I-TXK t Ml /ilX+2. 1 

J=TYKIO)/OY»2 .0 • 

KC< I . .1 1 =KG ( 1 . .11 + 1 

if(m.lt.io?oi fin m bi7 

MB I TF 1 14) (TXK(M) , TYK( N) . N= 1 , 1 020 ) 

M=Q 

817 K= K4- 1 

XP-XP-t-QVK 

816 COMTINJE 

Y SF P4P=Y5FPAR 6 1.0 

8010 KftAR=K 

KPA^KlAP/ 102-1 

MPPB = K-SAK-~- , PAR.-*10?C 

HPA«1=MPAR+1 J 

WRITE (14) t TXK< N ) »TYK (N ) tN-1 1 RPKR ) * 

' REWIND 14 

REWIND 15 


‘ C STORF PARTICLE CDORDINATFS ON TAPE 10 


IFJ_A8St TPP-T l.GT.ir.T/2.) 1 GC TO 834 

ASSIGN 834 TO KP c T 

0836 COKTMME 

PA TB( 1 ) =T 

DATBI? ) =K BAP 

DAT 6 ( 3 ) — I BAR 

DAT 8 ( 4 1 = J6AR 

DAT3I 5)=XR 

DAT6(5)=YT 

DATE ( 7 ) = 1 X 

OATq(?)*WY 

PATH 1)=1W 

DAT3(111=UVI 

DAT 8 ( 1 1 ) = VW 

DAT 5 [ 1 2 ) = PX 

DATtM I3)~=DY 

flATAI 14) = DTPP 

DAT fl< 151 = IWHEEL 

DATP(1S)=IFSYM 

DAT H ( 1 7 ) - F T AO 

0AT61 IB) -UWO 

GAT8( 191 = VHP 

DAT fU 7Q) = i-)*i 

OAT8I21 ) = 3miS 

1)ATB( ?2)-SHXSF 

DAT 6(23)-SDYSF 

Will IE ( 10) PATH 

MM =1071 

IIP 833 11 = 1 1 -3PAP 1 

IF( M .LE.MPAR ) OP Tfl 835 

MH = IIPPB 

0815 RT AD 111) ( IXMfl) i TYK1M) ,M=1,MM ) 

HR | IF 1 IOH TXMM). TYK(M) ,M=1,MM) 

083 A COMTIWJ c 

Fi'JD r ~I ■„ I; IP 

TPP=TPP+PTPP 

_NPP=FiPP+l 
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KiliODUCIBILITY OF THE 
ORIGINAL PA'-tt Tg P00R 



MRJ.T*- ((,. 26 1T.N PP 

RF W I NOJ 4 

GO T f 1~ <ii L T 

0034 com rp-iijg~ 

c 

c orr.iiiM oo- petcac cells 

c 

0820 no tin _j=2, jrari, 

on « i n l=~z% IBAR1 

LLLE i i . j i .nf.3 i r,u tii 050 ' 

I F 1 in C ( I.Jl.NE.O) GO TO (121 

F 1 1 t J ) =4 

phi ; -U.jj.f o. 

IF(n m.-lUNg.'i) GO Tu 830 

U.t..l.i.T ) ~0 .. _ . .. 

0630 ir<f( 1- 1 . J1.NS.4I Gfl Til 831 

111 I - 1 , 11 = 0. 

08 31 tF(F(l,.)+l).NE.4) GO TO R32 

V 1 1 .Jl = 0. 

0632 If lFH.J-H.Nc. 4) GO TCI 621 

VI 1 rJ-1 1-0. 

GO TO 621 

0350 1F1F1 l.J) .NF.-H GO Tfl 621 

PHI 1 t ,J)^0. 

IF1KC1 I.J 1 .F 0.01 on TO 821 

’ FI 1 , J 1 -3 

FlAlt . J ) = ET AO 

FC1 I . J 1 = 2 

0621 CONTINUE 

00 623 J=j,JBA»l 

DO 023 1=2.13431 

IF1KCI I, J ) .FO.O) GO TC 323 

IFI Ft 1, J1 .ME. 2 I GO Tfl 351 

IF! FI 1*1 . J l.Ng.4. AMD.ru- 1. J1.NE.4. AND.F1 1 .J + l 1 .NE.4.AN0. 

*FI 1. J-11.NF.41 Gil TO 6 23 

F( I.J ) = 3 

PHlIl.lbO, 

GO TO 823 

851 IF1FI I, J1 .IMF. 31 GO TO 823 . 

SlIPPHI^O. ^ 

N - Cl 

1F1F1 i + 1, J 1.E0.4.QR.FI I- l, J > .FQ.4.0R.F1 1 . J-t-1 ) ■ EC . 4. DR. 

«F1 I , J-ll.fcC.41 GO TCI 323 

1 F ( F I IH, Jl .ME. 2) GO TO 352 

SUNPHI = SUHPHH-PHI( t-Fl ,J) 

N = 1 n '<-1 

0852 1 F ( FI 1- 1 1 J ) .N£ . 2 I GO TC 353 

suNF u i=sunpmFpn[ t i-i.j) 

N=N + l 

0853 1F1F1 1.J + H.NE.2) GO TO 854 

SUHPHI=SUHPHI+PHI(1.J+11 

N-N 1 1 

0854 IF1FII.J— 11.NE.2) GO TO 855 . 

S1JMPHI=SUi3PHI»Pm ( 1 , J-l 1 

N=wn 

0855 Fll >01 =2 

FPN g N 

IF1N.EQ.0) GO TO 656 

PHI 1 I.Jl-SUNPHI/FPN 

GO Tl 623 

0856 IF! FE< I, J ) .EE .2) GO TP 623 

83ITE16,6)l,J,FU,J),fei 1, Jl.KClI.Jl 

0823 CONTINUE 


C Pc-lM?OS£ BOUNDAPY CONDITIONS TO FREE SURFACE 

C~ 

A5S1CN B«0 TC KRFT 

GO TO 545 

0800 CONTINUE 


C ESTABLISH ETA FOR POPPFRLY EMPTY CELLS 

C 

ASSIGN 824 TO KRFT 

no 024 j^a.jtiAPi 

DP 024 1=2,1 BAP L 

IF! FC1 I,J 1 .NC.2) 324,611 

0824 CONTINUE 

C " 

C TL-ST FOP TIME TP PRINT CELL VARIABLES 

C 

ASSIGN 801 TO KPET 

0802 COMINH? ' 

Mrfri 4 .ll) T.SFK..SFY 

WR1 T'r 16. IflPT , SDXSFf SDYSF 

Will if (6 . 22 I T. n-i sr, i:ysf 

W8 1 TF I . 30fT.UP~,V 3 , STNMAX 



HR I TF t6. 3 3IT.At,PHAU,ALPHAV,CnNHAX 

Wilts *•» * T . UN AX ,'l DM AXj j DMAX 

.i 2 LL'f. HH rm'l.eKKJRvrMJH 

TJPj.ff j 6.Ti) T.x sF {> Ap,y s fpah,k~bak 

hu'i rr <*,, fj i r, ic.CunV~,'rwk 

N; I li =j> 

tr'V A D < ;irCI > -T> . CT.(DT/?.l) go to SOI 

DA TAJ 1 1_^T 

DAfA("?~) -KOAR 

HAT A ( 3) =SOXSF 

OATA^MSOYSF 

DATA! 51-XSFPAR 

_D A TA m=YSf»AR 

PA T A ( 7) ~SFX 

DAT A 1 0 ) = S FY 

OA T A ( *J) - V H - 

DAT Al ni=i)W 

OA T A ( 11) =AI ^HAIJ 

D AT A( 12 I = Al PHAV 

DATA! 13) -Q TCP 

TlATAU-'. I^DTSP 

DATA 1 1t)=DTPP 

DAT A | Lftl^lftAP 

DATAt L?)=Ji>AR 

DATA(1»I1=XF 

data; ioi^vt 

PATA(20)=HX 

OATAI 21)=WY 

DAT A t 22 ) 

DATA ( 23 ) =DX 

DATA I 3A ) =DY 

PATA(25>=QXSF 

DATA I2M = PY SF 

OATAI27) = ERMnRij 

D4TA< ?d) ^EPRTiRV 

0 A T A ( 2v ) = I C 

DATA (30) = 5 DMA 

0ATA(31) = ’;FXD 

DATAI 32) = 9FYG 

DATAI33)=I^HFFL 

data; 3<ti = rFSv.-i 

OATA(33)=gTftO 

data; 3s>=uho 


DATA; 37) -VwO 

DATA ;38)=VP 

DATAI 39)=RH'DS 

HRITF(a)JATA»UiV 

HMD FILE 3 

WftI TF;i3l DATA,PHTt£TA,F.FE 

EMO FILE 13 

TCP^TCP+PTCP 

NCS=NGS+1 

• HRI FE (t> ,AQ30 ITtNCS 

GO TO <RE T 

0801 CONTE NUE 


C Efjn PROBLEM 


if; (TL-T1 .GT. (DT/2.1 1 GO TO 300 

0901 CALL SECOND ( fit IME) 

HR I TE ( a , 12 1 ST i ME 

GO TO 4000 

9909 CONTINUE 

STOP 

END 
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2 0(4.)),nuLS0Y(6O) ,H2 ( 60) , C ( ‘.0 ) 


fcl*SLN=L 
r r c i w) i 

PS 

5 *2,13 


2 N1»N-1 



O' I mi K=1 ,NCVS 


rr 


111 BHMBBBH 


UcLSOY ( I > = ( OFL-' ( 1 1 K ) -DPI Y ( ! IiK 1 1/ H2 III 


S2 I 1,K1=2.«»FLS0Y(I 


rMFEA- W071 7960 
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MlBIHTJin 


4 FF (T(.I)-XI 1)1 58, 17.55 
•55 IF (T[J)-X(N1 ) 57 , 59 , 5d 
56 IF 1TIJ1-XH )) 60.17.57 
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44_. FORMAT (I4,24HTH ARGUMENT OUT OF RANGE! 


lafliimaRwi 


HnwEnn 


oeta 


nufetHlIhrilHiflBMfai 


0ELS0S = (_S21 I i K J +■ S 2 t [I , K) +SS2 ( K , J ) ) / 5 


SStK.J )— Y t l.KH-HT K-OEL Y t I , K ) + PR0D<-QEL SOS 

SSI IK, J)=1)ELY( I.K l + IHTi + hT2}*l)ELSQSTPRni)*SX( I,KJ/6.0 


CKBllBiM 

n w aflgimnj i 


I F ( t rn.GT. 9 ) FFTUSi'l 


Oil 120 K- I I NC VS 


20 PROXIMK) -0.0 


62 PRCXINlK)=PRrXIN(K |t.5*HII) 


120 CONTIH'IE 


















TABLE I 

SUMMARY OF RESULTS 


Quantity 

Comp. 

'Exp. 

Comp. 

Exp. 


EXp. 

Equivalent mass of block, gm 

102.2 

102.2 

102.2 

102.2 

153 

153 

Initial impact velocity, cm/sec 

30 

33 

50 

47 

30 

30 

Peak acceleration, g 

1.76 

1.52 

2.84 

2.28 

1.24 

1.12 

Peak pressure, kN/m 2 

2.9 

6.3 

4.7 

9.2 

3.0 

6.5 

Maximum penetration, cm 

27.11 

23.07 

28.16 

24.46 

— 

— 

Time at maximum penetration, sec 

.395 

.390 

.375 

.373 

— 
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Strain-rate, 


Figure 1. ' Soil constitutive 
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•REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 



Figure 6. Experimental apparatus used to measure 

displacement, acceleration, and pressure. 





c. Figure 7. Impact block. 
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Pressure, kN/m 2 Acceleration, g- units Displacement, cm 




(a) v q ss 30 cm/ sec, mass of block ~ 102.2 gm 


Figure g. Comparison of computational and 
experimental time histories of 
rigid block impacting on still 
water . 


(a) Displacement of Block 


OOOO 



.1 .2 .3 .4 

(b) Deceleration of Block 



.1 .2 .3 .4 


Time, seconds 

'50 cm/sec, mass' of' block = 102, 2-sgm 


Figure ' 8 . Continued . 




(a) Displacement of Block 
to 



(b) Deceleration of Block 

<N 



Time, seconds 


(c) v q = 30 cm/sec, mass of block = 153 gm 
Figure 8. Concluded. 


70 


'O ^ 

|8 

1 | 

l| 

>— i ►<! 

M 

TJ 

o !* 

© f-9E 

» W 
KS 


10 


20 


31 


41 


62 



10 


Field -width. cm. 
20 31 41 


62 



(b) Particle configuration, t = .155 sec. 


16 


30 


45 



(c) Velocity vector configuration^ = .1S5 sec. 


Computational Parameters 


Ax =Ay = .9524cm. 
No. of cells =62x62 
v o * 50 cm/sec, 

Mass =102. 2gm. 


Initial fluid depth = 24 cells (23cm) 

Width of fluid (half field) =60cells (57cm) 

Size of half block =4 x 34 cells! 3. 8 x 32.4cm) 
Fluid kinematic viscosity = .01 cm 2 / sec. 


Figure 9. Fluid dynamics evolution for typical 
computational study. 
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Pressure, kN/m 2 


Figure 10. Computed excess hydrostatic pressure at center of 
base of rigid block .at' time t = .003 seconds 
after impact. Mass — 102.2 gm 
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Figure 11. Time, history of forces on 
rigid block of various 
•contact length. 
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Figure 12 . Effect of contact len< 
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Figure 13. Effect of initial horizontal block velocity 
on forces. 
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